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;

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

[Krjuger]

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

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

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

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

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

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

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

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