Автор: Maximys 1.5.2010, 7:03
У меня есть задача по Численным методам: написать программу для решения уравнений методом неполной релаксации, наискорейшего спуска и минимальных невязок. Алгоритм по методу неполной релаксации успешно реализован, но два других метода - просто "труба". Нашел в сети методичку СамарскогоГосУнивера, там для метода наискорейшего спуска даются следующие формулы: 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 в формулах означает транспонирование матрицы, стоящей слева от него