Проблема с реализацией метода наискорейшего спуска для решения СЛАУ

44K

26 августа 2010 года

18 / / 01.01.2010

Необходимо на языке ObjectPascal реализовать метод Наискорейшего спуска для решения Систем Линейных Алгебраических Уравнений(СЛАУ). Нашел необходимый материал(книги под авторством Фадеевых; Петрова, Лобанова и т.д.), применил формулы, указанные там и вроди бы при "хороших" СЛАУ все отлично работает - система, решенная в книге Фадеевых решается точно также(одинаковые невязки, корни), некоторые другие системы тоже отлично решаются, но когда пытаюсь решить СЛАУ со следующими матрицами коэффициентов и свободных членов, корни "улетают" в бесконечность. Проблемная матрица A и вектор свободных членов B имеют следующий вид:
(1.7 2.8 1.9)
A=(2.1 3.4 1.8)
(4.2 -1.7 1.3)

(0.7)
B=(1.1)
(2.8)
Как известно, метод наискорейшего спуска применим для любой системы, матрица коэффициентов которой имеет определитель больше нуля. Переставляем строки матрицы A, чтобы на главной диагонали стояли максимальные по модулю элементы. Делим коэффициенты и свободные члены каждого уравнения на диагональные элементы, в результате получается положительно-определенная матрица - но почему при её решении получаются заоблачные корни и программа "вылетает":confused:? Подскажите, пожалуйста, кто что сможет по этому поводу:(.
Код моей процедуры следующий:

Свернуть исходник

Код:

procedure MethodOfSteepestGradientDescent(U:SLAU; var Ot:ROOTS; aSize:integer; NeedAccuracy:real; var ColItterationG:longint);

type ArrayWithNevazGD = array[1..MaxSize] of real;

  const zero : integer = 0;

  const one : integer = 1;

    var i : integer;

       li, lj:integer;

        CurCalc : real;

        R : ArrayWithNevazGD;

        UMulR : ArrayWithNevazGD;

        m : real;

        MaxNev : real;

   //локальная функция для перемножения массивов типа ArrayWithNevaz

  function MulVectorToVector(Any1, Any2: ArrayWithNevazGD; aaSize:integer) : real;

      var lli : integer;

          OutDigit : real;

  begin

    OutDigit:=zero;

    for lli:=one to aaSize do

      OutDigit:=OutDigit + Any1[lli]*Any2[lli];

  MulVectorToVector:=OutDigit;

  end;

  //---------------------------------------------------------------

  //локальная функция для поиска максимальной невязки в массиве типа ArrayWithNevaz

  function FindMaxNevaz(aR : ArrayWithNevazGD; aaSize:integer) : real;

      var lllli : integer;

          OutMaxNevaz: real;

  begin

  OutMaxNevaz:=zero;

  for lllli:=one to aaSize do

    if(abs(aR[lllli])>OutMaxNevaz) then

        OutMaxNevaz:=abs(Ar[lllli]);

  FindMaxNevaz:=OutMaxNevaz;

  end;

  //-------------------------------------------------------------------------------

begin

{формула для метода Наискорейшего спуска:

    Ot[p+1] = Ot - m*R

      m = (R, R) / (A*R, R)}

ColItterationG:=zero;

for li:=one to aSize do begin

      CurCalc:=zero;

      for lj:=one to aSize do

        CurCalc:=CurCalc + U[li,lj] * Ot[lj];

      R[li]:=U[li, aSize + one] - CurCalc;

    end;

repeat

  inc(ColItterationG);

for li:=one to aSize do begin

     UMulR[li]:=zero;

     for lj:=one to aSize do

       UMulR[li]:=UMulR[li] + U[li, lj] * R[lj];

    end;

m:=MulVectorToVector(R, R, aSize) / MulVectorToVector(UMulR, R, aSize);

  for i:=one to aSize do

     Ot:=Ot + m*R;

for li:=one to aSize do begin

      CurCalc:=zero;

      for lj:=one to aSize do

        CurCalc:=CurCalc + U[li,lj] * Ot[lj];

      R[li]:=U[li, aSize + one] - CurCalc;

    end;

  MaxNev:=FindMaxNevaz(R, aSize);

until MaxNev<NeedAccuracy;

end;

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