Математика \ Методы вычислений

Минимизация квадратных функций. Градиентный метод наискорейшего спуска, страница 2

где r(k)=Ax(k)-b, αk=(r(k),r(k))/(Ar(k),r(k))

Как мы увидим, при сделанных предположениях о матрице А последовательность {x(k)} сходится к решению системы (2).

Рассмотрим количество операций умножения, деления и сложения: для r(k)=Ax(k)-b –2n^2+ n опер.;  для αk=(r(k),r(k))/(Ar(k),r(k)) – 4n +2n^2 опер.;

для x(k)= x(k)-α0r(k) – 2n опер. Таким образом для каждого шага 7n+4n^2 опер.

ТЕОРЕМА 1Последовательность приближения метода наискорейшего спуска {x(k)} для системы (2) сходятся к её решению, при этом имеет место оценка

|| x(m) - x*|| =<(||r(0)||/ λn)[( λ1-λn)/( λ1+ λn)]m , m=1,2…

где λ1-наибольшее и λn – наименьшее собственные числа матрицы А и r(0)=A x(0)-b.

Таким образом из этой теоремы мы можем приблизительно сказать на каком шаге (m) остановится наш процесс приближения, для этого надо рассмотреть уравнение

(||r(0)||/λn)[(λ1-λn)/(λ1+λn)]m=0, с неизвестной m, которая напрямую зависит от начального приближения.

ОТВЕТ задание 2

Программа написана на delphi5.0

КОД:

НЕОБХОДИМЫЕ ФУНКЦИИ:

//функция преобразование строки в целочисленный тип

function si(t:string):integer;

begin

si:=strtoint(t);

end;

//функция преобразование строки в вещественный тип

function sf(t:string):real;

begin

sf:=strtofloat(t);

end;

//функция преобразование вещественного числа в строковый тип

function fs(t:real):string;

begin

fs:=floattostr(t);

end;

//функция возведения в квадрат

function kv(t:real):real;

begin

kv:=t*t;

end;

ПЕРЕМЕННЫЕ:

//двумерный массив размерности 3600(матрица А), вещественный тип

a: array [1..60,1..60] of real;

//члены последовательности {xk}

x: array [-1..10000,1..4] of real;

//члены последовательности {rk}

r: array [0..10000,1..4] of real;

//одномерные массивы

b,ax,ar,axb,xx: array [1..4] of real;

//переменные - вещественный тип

z,z1,f1,f2: real;

//переменные - целочисленный тип

i,j,k,d:integer;

f: TextFile;//переменная - тип текстового файла

fName: String[80]; //переменная - для названия файла

buf: String[80]; //переменная - буфер для хранения текста из файла

text1:string; // "кусочек" буфера

ПРОЦЕДУРА АВТОЗАПОЛНЕНИЯ ДАННЫХ НАХОЖДЕНИЯ МИНИМУМА:

//автозаполнение из файла

if radiogroup1.ItemIndex=0 then begin

fName := 't16-w.txt'; //название файла

AssignFile(f, fName);//открытие этого файла

{$I-}   //директори для проверки правильности открытия

Reset(f);

{$I+}

//сама проверка

if IOResult <> 0 then

begin

MessageDlg('Ошибка доступа к файлу ' + fName,

mtError,[mbOk],0);

exit;

end;

k:=0;

text1:='';

while not EOF(f) do //цикл для считывания строк файла

begin

k:=k+1;

d:=0;

readln(f, buf); //считывание 1-й строки

for i:=1 to length(buf) do begin //цикл для заполнения массива

if (copy(buf,i,1)=' ')  then else begin

if copy(buf,i,1)='.' then text1:=text1+',' else

text1:=text1+copy(buf,i,1);

end;

if (copy(buf,i,2)='  ') or (i=length(buf)) then begin

a[d+1][k]:= sf(text1);

text1:='';

d:=d+1;

end;

end;

end;

CloseFile(f); //закрытие файла

st2.RowCount:=d-1; //таблица начального приближения размерности oxn

st3.RowCount:=d-1; //таблица ответа

st4.RowCount:=d-1;

end;

if radiogroup1.ItemIndex=1 then begin

//нахождение минимума функционала

k:=0;

d:=st3.RowCount; //размерность исходных данных

//считывание данных в массивы

for i:=1 to d do begin

b[i]:= -a[d+1][i]; //вектор b

x[0,i]:= sf(st2.Cols[0][i-1]); //вектор х(начальное приближение)

x[-1,i]:= 10;

end;

f1:=2;

f2:=1;

// xx[d] - ||xk-xk-1||^2 (один из критериев остановки спуска)

xx[d]:=kv(x[k,1]-x[k-1,1])+kv(x[k,2]-x[k-1,2])+kv(x[k,3]-x[k-1,3])+kv(x[k,4]-x[k-1,4]);

// axb[d] - ||f'(xk) ||^2 (один из критериев остановки спуска)

axb[d]:=1;

while (sqrt(xx[d])>exp(ln(10)*(-15))) //второй критерий

and (abs(f1-f2)>exp(ln(10)*(-15))) //третий критерий

and (sqrt(axb[d])>exp(ln(10)*(-15))) {первый критерий} do

begin

//цикл нахождения А*хk

for j:=1 to d do begin

ax[j]:=0;

for i:=1 to d do  ax[j]:= ax[j]+a[i,j]*x[k,i];

end;

//цикл нахождения rk

for i:=1 to d do  r[k,i]:= ax[i]- b[i];

//цикл нахождения A*rk

for j:=1 to d do begin

ar[j]:=0;

for i:=1 to d do  ar[j]:= ar[j]+a[i,j]*r[k,i];

end;

z:=0;

z1:=0;

//цикл нахождения ak=z/z1

for i:=1 to d do

begin

z:=z+r[k,i]*r[k,i];

z1:=z1+ ar[i]*r[k,i];

end;

//цикл нахождения xk+1

for i:=1 to d do x[k+1,i]:=x[k,i]-(z/z1)*r[k,i];

k:=k+1;

//цикл нахождения значения функционала в т. xk-1

if k>3 then begin

z:=0;

z1:=0;

Скачать файл