LU разложение. Опции

[Krjuger]

Собственно есть задача, полную ее суть писать не стану,но проблема заключается в следующем нужно с помощью LU разложения решить систему.Рабочий код на С++ у меня был и я его попытался перенести на Delphi.Но к несчастью моя затея обернулась неудачей.

type
  TVec = array of Extended; // векторный тип
  TMatr = array of TVec; // матричный тип

function  LUrazl(matr: Tmatr): Tvec;
var
  i, j,k, row, col: Word;
  res: TVec;
  _matr_: TMatr;
  L,U: TMatr;
  Sum: double;
begin
  _matr_ := CopyMatr(matr); // Copy не работает как надо
	row := Length(_matr_);
  col := Length(_matr_[0]);
  SetLength(L,row,col);
  SetLength(U,row,col);
	if row <> col-1 then
    ShowMessage('Ошибка: некорректные размерности матрицы! (МГ)'); //error
  for i:=0 to  row-1 do
    for j:= 0  to col-1 do
   begin
      L [i][j] := 0;
      U [i][j] := 0;

      if i = j then
        L [i][j] := 1;
   end;



//==============================================

//находим первый столбец L[][] и первую строку U[][]

  for i:= 0 to row-1 do
  begin
    L [i][0]:= _matr_ [i][0]/ _matr_ [0][0];
    U [0][i]:= _matr_ [0][i] / L [0][0];
  end;

    for i:= 0 to row-1 do
     begin
          for j:= 0 to col-1 do
            begin
                   U [0][ i] := _matr_[0][ i];
                   L [i][ 0] := _matr_[i][ 0] / U[0][ 0];
                   sum := 0;
                   for k:= 0 to i  do
                   begin
                       sum :=sum+ L[i][ k] * U[k][ j];
                   end;
                   U[i][ j] := _matr_[i][ j] - sum;
                   if (i > j) then
                     L[j][ i] := 0
                   else
                    begin
                       sum := 0;
                       for  k:= 0 to i do
                       begin
                           sum:=sum+L[j][ k] * U[k][ i];
                       end;
                       L[j][ i] := (_matr_[j][ i] - sum) / U[i][ i];//В этой строке выдает деление на ноль так же в дебаге видно,что значение i=65000+
                    end;

            end;
          res[i]:= U[i][ i];
     end;
end;

Непонятна причина такого поведения,ведь на С++ все это работало.Если надо могу добавить файл целиком.

[IUnknown]

Добавь лучше С++-ную реализацию. Посмотрим, что сделано не так, почему не работает...

[Krjuger]

Вот С++,так как оно мне требовалось, чтобы посчитать стопочку типовых расчетов,особого энтузиазма в оформлении нет.

#include <iostream>
using namespace std;

int main ()
{
  const int n=4;
  double sum = 0;
  double A [ n ][ n ]= { 
{ 2, -5,6,1},
{ -16, 48,-55,-5},
{ 0, 24,-25,10},
{8,21,30,4 }};
  double L [ n ][ n ];
  double U [ n ][ n ];

//задаем матрицу A[][] ...

  for (int i  = 0; i < n; i++)
  {
   for (int j = 0; j < n; j++)
   {
 //     cout << "\na[" << i+1 << "][" << j+1 << "] = ";
   //   cin >> A [i][j];
     L [i][j] = 0;
     U [i][j] = 0;

      if (i == j)
        L [i][j] = 1;
   }
  }

//==============================================

//находим первый столбец L[][] и первую строку U[][]

  for (int i = 0; i < n; i++)
  {
    L [i][0] = A [i][0]/ A [0][0];
    U [0][i] = A [0][i] / L [0][0];
  }

for (int i = 0; i < n; i++)
           {
               for (int j = 0; j < n; j++)
               {
                   U [0][ i] = A[0][ i];
                   L [i][ 0] = A[i][ 0] / U[0][ 0];
                   double sum = 0;
                   for (int k = 0; k < i; k++)
                   {
                       sum += L[i][ k] * U[k][ j];
                   }
                   U[i][ j] = A[i][ j] - sum;
                   if (i > j)
                   {
                       L[j][ i] = 0;
                   }
                   else
                   {
                       sum = 0;
                       for (int k = 0; k < i; k++)
                       {
                           sum += L[j][ k] * U[k][ i];
                       }
                       L[j][ i] = (A[j][ i] - sum) / U[i][ i];
                   }
               }
           }


//====================================================
   cout << "\n\n";
   for (int i = 0; i < n; i++)
   {
     for (int j = 0; j < n; j++)
	 {	
		 if (i == j)
          L [i][j] = 1;
       cout << "  " << L [i][j] << "   ";
	 }
     cout << "\n\n";
   }
   

   cout << "\n\n";

   for (int i  = 0; i < n; i++)
   {
     for (int j = 0; j < n; j++)
       cout << "  " << U [i][j] << "   ";
     cout << "\n\n";
   }

   return 0;
}

Так как типовые были сдаты, и пара вариантов были посчитаны руками и результаты были одинаковыми,поэтому я счел,что он работает верно.

Да и забыл упомянуть для текущей задачи мне требуеться,чтобы оно считалось для максимальных размерностей 101х100.А то у меня уже был случай,когда Гаусс верно считал для любой матрицы меньше 20,а для больших переставал работать,пришлось полностью пересмотреть концепцию)))

Сообщение отредактировано: Krjuger - 1.11.2011 18:18

[IUnknown]

У тебя как минимум в двух местах неправильно сделан перенос с С++... Оба - циклы по K. И там и там надо ходить до (i - 1), а не до i... Тестировал программу на той же матрице, что и Сишный вариант?

P.S.

Матрицы вот такие должны получаться (это результат работы твоего чуть-чуть подпиленного дельфийского кода):

L-matrix:
    1.0000    0.0000    0.0000    0.0000
   -8.0000    1.0000    0.0000    0.0000
    0.0000    3.0000    1.0000    0.0000
    4.0000    5.1250  -10.4688    1.0000
U-matrix:
    2.0000   -5.0000    6.0000    1.0000
    0.0000    8.0000   -7.0000    3.0000
    0.0000    0.0000   -4.0000    1.0000
    0.0000    0.0000    0.0000   -4.9063

?

Да, кстати, у тебя там балаган с возвращаемыми значениями, ты под Res память-то не выделил. Да и как собрался возвращать и L- и U-матрицы? Проще передать через параметры...

Сообщение отредактировано: IUnknown - 1.11.2011 19:20

[Krjuger]

Нет не тестировал,я смысла не вижу,если оно даже не отрабатывает до конца,то говорить правильности решения пока рано.

У тебя как минимум в двух местах неправильно сделан перенос с С++... Оба - циклы по K. И там и там надо ходить до (i - 1), а не до i...

К несчастью так нельзя,когда мы самый первый раз заходим в цикл по i у нас i=0, и i-1 в цикле по K мне даст большую беду,в Сях то такой цикл просто проскочился бы не отробатывая,а тут так не катит(((

Да, кстати, у тебя там балаган с возвращаемыми значениями, ты под Res память-то не выделил. Да и как собрался возвращать и L- и U-матрицы? Проще передать через параметры...

Мне не нужно возвращать эти матрицы,мне нужно вернуть вектор, который состоит из элементов диагонали матрици U это и будут решения системы,такой уж метод.

Матрицы вот такие должны получаться (это результат работы твоего чуть-чуть подпиленного дельфийского кода):

В Сях то он таким и получаеться.

Сообщение отредактировано: Krjuger - 1.11.2011 19:26

[IUnknown]

К несчастью так нельзя,когда мы самый первый раз заходим в цикл по i у нас i=0, и i-1 в цикле по K мне даст большую беду

Бла-бла-бла... Кроме Word-а в Дельфи есть еще и Integer, который способен принимать отрицательные значения. Чуешь, откуда ошибка?

В общем, вот, прогнал на FPC твой дельфийский код, ничего не вылетает, печатает то, что приведено выше:

type
   TVec = array of Extended;
   TMatr = array of TVec;

procedure LUrazl(_matr_: Tmatr; var L, U : TMatr);
var
   i, j, row, col: Word;
   k : integer;
   res: TVec;
   Sum: double;
begin
   // _matr_ := CopyMatr(matr);
	row := Length(_matr_);
   col := Length(_matr_[0]);
   SetLength(L,row,col);
   SetLength(U,row,col);
   if row <> col-1 then
      writeln('error');

   for i:=0 to  row-1 do
      for j:= 0  to col-1 do
      begin
         L [i][j] := 0;
         U [i][j] := 0;
         if i = j then L [i][j] := 1;
      end;

   //==============================================

   for i:= 0 to row-1 do
   begin
      L [i][0]:= _matr_ [i][0]/ _matr_ [0][0];
      U [0][i]:= _matr_ [0][i] / L [0][0];
   end;

   for i:= 0 to row-1 do
   begin
      for j:= 0 to col-1 do
      begin
         U [0][ i] := _matr_[0][ i];
         L [i][ 0] := _matr_[i][ 0] / U[0][ 0];
         sum := 0;
         for k:= 0 to i-1  do
         begin
            sum :=sum+ L[i][ k] * U[k][ j];
         end;
         U[i][ j] := _matr_[i][ j] - sum;
         if (i > j) then
            L[j][ i] := 0
         else
         begin
            sum := 0;
            for  k:= 0 to i-1 do
            begin
               sum:=sum+L[j][ k] * U[k][ i];
            end;
            L[j][ i] := (_matr_[j][ i] - sum) / U[i][ i];
         end;
      end;
   end;
end;

var
  A : tmatr;
  L, U : tMatr;
  i, j : integer;

begin
  setlength(a, 4, 4);

  a[0, 0] :=   2; a[0, 1] := -5; a[0, 2] :=   6; a[0, 3] :=  1;
  a[1, 0] := -16; a[1, 1] := 48; a[1, 2] := -55; a[1, 3] := -5;
  a[2, 0] :=   0; a[2, 1] := 24; a[2, 2] := -25; a[2, 3] := 10;
  a[3, 0] :=   8; a[3, 1] := 21; a[3, 2] :=  30; a[3, 3] :=  4;

  LURazl(A, L, U);

  writeln('L-matrix:');
  for i := 0 to Length(L) - 1 do
  begin
     for j := 0 to Length(L[0]) - 1 do
        write(L[i, j]:10 :4);
     writeln;
  end;
  writeln('U-matrix:');
  for i := 0 to Length(U) - 1 do
  begin
     for j := 0 to Length(U[0]) - 1 do
        write(U[i, j]:10 :4);
     writeln;
  end;

end.

Дельфи запускать лень

Сообщение отредактировано: IUnknown - 1.11.2011 19:26

[-Татьяна-]

Собственно есть задача, полную ее суть писать не стану,но проблема заключается в следующем нужно с помощью LU разложения решить систему.Рабочий код на С++ у меня был и я его попытался перенести на Delphi.Но к несчастью моя затея обернулась неудачей.

type
  TVec = array of Extended; // векторный тип
  TMatr = array of TVec; // матричный тип

function  LUrazl(matr: Tmatr): Tvec;
var
  i, j,k, row, col: Word;
  res: TVec;
  _matr_: TMatr;
  L,U: TMatr;
  Sum: double;
begin
  _matr_ := CopyMatr(matr); // Copy не работает как надо
	row := Length(_matr_);
  col := Length(_matr_[0]);
  SetLength(L,row,col);
  SetLength(U,row,col);
	if row <> col-1 then
    ShowMessage('Ошибка: некорректные размерности матрицы! (МГ)'); //error
  for i:=0 to  row-1 do
    for j:= 0  to col-1 do
   begin
      L [i][j] := 0;
      U [i][j] := 0;

      if i = j then
        L [i][j] := 1;
   end;
//==============================================

//находим первый столбец L[][] и первую строку U[][]

  for i:= 0 to row-1 do
  begin
    L [i][0]:= _matr_ [i][0]/ _matr_ [0][0];
    U [0][i]:= _matr_ [0][i] / L [0][0];
  end;

    for i:= 0 to row-1 do
     begin
          for j:= 0 to col-1 do
            begin
                   U [0][ i] := _matr_[0][ i];
                   L [i][ 0] := _matr_[i][ 0] / U[0][ 0];
                   sum := 0;
                   for k:= 0 to i  do
                   begin
                       sum :=sum+ L[i][ k] * U[k][ j];
                   end;
                   U[i][ j] := _matr_[i][ j] - sum;
                   if (i > j) then
                     L[j][ i] := 0
                   else
                    begin
                       sum := 0;
                       for  k:= 0 to i do
                       begin
                           sum:=sum+L[j][ k] * U[k][ i];
                       end;
                       L[j][ i] := (_matr_[j][ i] - sum) / U[i][ i];//В этой строке выдает деление на ноль так же в дебаге видно,что значение i=65000+
                    end;

            end;
          res[i]:= U[i][ i];
     end;
end;

Непонятна причина такого поведения,ведь на С++ все это работало.Если надо могу добавить файл целиком.

А ВЫ не могли бы сбросить вашу задачу пожалуйста, она мне очень нужна...(lovetanya.15@mail.ru) Спасибо