У меня есть задача по Численным методам: написать программу для решения уравнений методом неполной релаксации, наискорейшего спуска и минимальных невязок. Алгоритм по методу неполной релаксации успешно реализован, но два других метода - просто "труба". Нашел в сети методичку СамарскогоГосУнивера, там для метода наискорейшего спуска даются следующие формулы: X[k+1] = X[k] - a[k] * (b - A*X[k]), где параметр a[k] ищется по формуле: a[k]=((A*X[k] - B )T*(A*X[k] - B )) / ((A*X[k] - B )T*A*(A*X[k] - B )T). Применяю эти формулы в своей программе - метод сперва сходится(приближается к корням СЛАУ), а затем "проскакивает" и корни начинают "уходить в бесконечность"(метод проходит "мимо" корней СЛАУ). К примеру, начальные приближения корней равны x1=1; x2=3; x3=2.3; корни СЛАУ равняются x1=0.5; x2=1.34; x3=2; метод доходит до значения корней x1=0.75; x2=2.33; x3=1.8(т.е. приближается к действительным корням СЛАУ), а затем корни равны x1=2; x2=3.4; x3=3...("постепенно уходим далеко от корней"). У меня имеются большие сомнения насчет корректности формул из методички - возможно к корням надо приближаться последовательно(увеличивать/уменьшать сперва один корень, затем второй и т.д.)? Вот мой программный код(синтаксис языка ObjectPascal с использованием локальных процедур и функций):
Код

procedure MethodOfSteepestGradientDescent(U:Uravnenia; var O:Otvet; KU:integer; NeedAccuracy:real);
type VectorNevazok=array[1..10] of real; //локальный тип данных для вектора невязок
    var i,j,k:integer;
        R:VectorNevazok;
        ResultMulVectorToMatrix:VectorNevazok;
        numerator:real;
        denominator:real;
        Alpha:real;
//локальная функция для подсчета невязок на k-ом шаге
function Discrepancies(Ur:Uravnenia; AnswerToStep:Otvet; NumberOfEquations:integer):VectorNevazok;
    var iL,jL:integer;
        NowDiscrepancies:real;
begin
for iL:=1 to NumberOfEquations do begin
     NowDiscrepancies:=0;
     for jL:=1 to NumberOfEquations do
         NowDiscrepancies:=NowDiscrepancies + Ur[iL,jL]*AnswerToStep[1, jL]; //r=a[i,1]*x1 + a[i,2]*x2 + ... + a[i,NumberOfEquations]
     NowDiscrepancies:=NowDiscrepancies - Ur[iL, NumberOfEquations + 1];//r=(левая часть уравнения) - (правая часть уравнения)
     Discrepancies[iL]:=NowDiscrepancies;
   end;
end;
//==============================================

//локальная процедура для сравнения вектора невязок с допустимой невязкой
  function CanProcedureStop(Nevaz:VectorNevazok; HowManyEquations:integer; NeedTochnost:real):boolean;
      var iC:integer;
          LocalBool:boolean;
  begin
  LocalBool:=true;
  iC:=1;
  while(iC<=HowManyEquations)and(LocalBool=true) do begin
      if(abs(Nevaz[iC])>NeedTochnost) then
         LocalBool:=false;
      inc(iC);
    end;
  CanProcedureStop:=LocalBool;
  end;
  //=============================================

//локальная функция для умножения вектора Vec на матрицу Mat1
function MulVectorToMatrix(Vec:VectorNevazok; Mat1:Uravnenia; m:integer):VectorNevazok;
    var iM,jN,kM:integer;
        rMul:real;
begin
for iM:=1 to m do
   for jN:=1 to m do begin
      rMul:=0;
      for kM:=1 to m do
        rMul:=rMul + Mat1[kM, iM]*Vec[kM];
      MulVectorToMatrix[iM]:=rMul;
    end;
end;
//==================================================

//----локальная функция для умножения вектора-строки на вектор-столбец с размерностями w
function MulVectorLineToVectorColumn(VectorLine:VectorNevazok; VectorColumn:VectorNevazok; w:integer):real;[/b]
    var kMV:integer;
        rMV:real;
begin
rMV:=0;
for kMV:=1 to w do
     rMV:=rMV + VectorLine[kMV]*VectorColumn[kMV];
MulVectorLineToVectorColumn:=rMV;
end;
// ================================================================================
============

//сама процедура для уточнения корней системы уравнений
procedure NewRoot(QuantityUrav:integer; alpha:real; Nevaz:VectorNevazok; var Roots:Otvet);
      var iNR, jNR:integer;
  begin
  for iNR:=1 to QuantityUrav do
      Roots[1, iNR]:=Roots[1, iNR] - alpha*Nevaz[iNR];
end;
//==========================================================

begin
//Основная формула метода наискорейшего градиентного спуска:
//x[k+1]=x[k] - Alpha[k]*gradF(x[k]) = x[k] - Alpha[k]*(b - A*x[k]) = x[k] + Alpha[k]*(A*x[k] - b)
//в которой Alpha[k] ищется как:
//Alpha[k]=(((A*x[k] - b)T) * (A*x[k] - b))/((A*x[k] - b)T) * A * (Ax[k] - b))
repeat
  R:=Discrepancies(U, O, KU);
  numerator:=MulVectorLineToVectorColumn(R, R, KU); //т.к. вектор R имеет KU строк и один столбец(он является одномерным массивом), то транспонировать его не надо
  ResultMulVectorToMatrix:=MulVectorToMatrix(R, U, KU);
    denominator:=MulVectorLineToVectorColumn(ResultMulVectorToMatrix, R, KU);
  Alpha:=numerator/denominator;
  NewRoot(KU, Alpha, R, O);
until CanProcedureStop(R, KU, NeedAccuracy)



end;

Подскажите пожалуйста, что необходимо исправить в моём коде или хотя бы объясните, в чем ошибка?
В книге Петрова, Лобанова - "Лекции по Вычислительной математике" нашел информацию о методе минимальных невязок. Формулы для него следующие: X[k+1]=X[k] - t[k]*r[k], где t[k]=(A*r[k],r[k])/(A*r[k], A*r[k]). Если правильно понял, то (A*r[k],r[k]) означает векторное произведение, но где Петров и Лобанов видели формулы, к примеру, для векторного произведения 3-порядка(пространства), 4-порядка и т.д.? Поправте меня пожалуйста, если я что-то путаю в выражениях (A*r[k],r[k]) и (A*r[k], A*r[k])(если это не векторное произведение) и подскажите пожалуйста, чем тогда это является.
Символ T в формулах означает транспонирование матрицы, стоящей слева от него