Метод Дэвидона - Флетчера - Пауэлла

Первоначально метод был предложен Дэвидоном и затем развит Флетчером и Пауэллом. Метод Дэвидона-Флетчера-Пауэлла называют также и методом переменной метрики. Он попадает в общий класс квазиньютоновских процедур, в которых направления поиска задаются в виде -Dj*grad(f(y)). Направление градиента является, таким образом, отклоненным в результате умножения на -Dj, где Dj - положительно определенная симметрическая матрица порядка nxn, аппроксимирующая обратную матрицу Гессе. На следующем шаге матрица Dj+1 представляется в виде суммы Dj и двух симметрических матриц ранга один каждая. В связи с этим схема иногда называется схемой коррекции ранга два.

Фрагментыисходногокода.

#include"stdafx.h"

#include<math.h>

#include<iostream>

#include<conio.h>

usingnamespace std;

typedefdouble coord[2];

double fn(coord x)

{

return 2*x[0]*x[0]+4*x[0]*x[1]+x[1]*x[1];

}

double f1(double _lam, coord _x, coord _p)

{

coord x;

coord x1,p1;

p1[0] = _p[0];

p1[1] = _p[1];

x1[0] = _x[0];

x1[1] = _x[1];

x[0] = x1[0]+_lam*p1[0];

x[1] = x1[1]+_lam*p1[1];

return fn(x);

};

double goldenSection(coord _x, coord _p, double _eps)

{

constdouble c1 = (3-sqrt(5.0))/2;

constdouble c2 = (sqrt(5.0)-1)/2;

coord x;

x[0]=_x[0];

x[1]=_x[1];

double a = 0;

double b = 2;

double x1, x2;

x1 = a + c1*(b-a);

x2 = a + c2*(b-a);

while (c2*abs(b-a) > _eps)

{

if (f1(x1,x,_p) > f1(x2,x,_p))

{

a = x1;

x1 = x2;

x2 = a + c2*(b-a);

}

else

{

b = x2;

x2 = x1;

x1 = a + c1*(b-a);

}

}

return (a+b)/2;

};

double dfdx(coord _x, int _arg)

{

if (_arg == 0)

return 4*_x[0]+_x[1];

else

return _x[1]+2*_x[0];

}

int _tmain(int argc, _TCHAR* argv[])

{

setlocale(LC_ALL,"Russian");

double eps;

coord x;

coord res;

//cout<<"ВведитеХ0 "<<'\n'; cin>>x[0]>>x[1];

//cout<<"Введите точность"<<'\n'; cin>>eps;

x[0]=-230;

x[1]=-230;

eps=0.01;

int it = 0;

double H[2][2] =

{ {4, 1}, //Матрица H - единичная

{1, 2} };

double A[2][2], B[2][2], UUT[2][2], tempB[2][2];

double alf;

coord p,gradf,gradf0,U,V,tempHU;

gradf0[0] = dfdx(x,0); //Вычислили градиент в т. X

gradf0[1] = dfdx(x,1);

p[0] = - (H[0][0]*gradf0[0]+H[0][1]*gradf0[1]); //Нашли новое p = -H*f'(x)

p[1] = - (H[1][0]*gradf0[0]+H[1][1]*gradf0[1]);

while (sqrt(p[0]*p[0]+p[1]*p[1])>eps)

{

it++;

alf = goldenSection(x,p,1e-3); //Вычислилиальфа

x[0] = x[0] + alf*p[0];

x[1] = x[1] + alf*p[1]; //Нашли новое X

gradf[0] = dfdx(x,0); //Значение в градента в найденной точке

gradf[1] = dfdx(x,1);

V[0] = alf*p[0]; //V = Xi+1- Xi

V[1] = alf*p[1];

U[0] = gradf[0]-gradf0[0]; //U = f'(Xi+1)- f'(Xi)

U[1] = gradf[1]-gradf0[1];

double VU = V[0]*U[0]+V[1]*U[1];

for (int i=0; i<2; i++)

{

for (int j=0; j<2; j++)

{

A[i][j] = V[i]*V[j]/VU; //A = V*Vt/<U,V>

UUT[i][j] = U[i]*U[j]; //Вычисляемматрицу U*Ut

}

}

tempHU[0] = H[0][0]*U[0]+H[0][1]*U[1];

tempHU[1] = H[1][0]*U[0]+H[1][1]*U[1]; //Dычислили H*U

double HUU = tempHU[0]*U[0]+tempHU[1]*U[1]; //Вычислили<H*U,U>

for (int i =0; i<2; i++)

{

for (int j =0; j<2; j++)

{

tempB[i][j] = 0;

for (int r=0; r<2; r++)

{

tempB[i][j] += H[i][r]*UUT[r][j]; //tempB = H*U*Ut

}

}

}

for (int i =0; i<2; i++)

{

for (int j =0; j<2; j++)

{

B[i][j] = 0;

for (int r=0; r<2; r++)

{

B[i][j] += tempB[i][r]*H[j][r]/HUU; //B = H*U*Ut*Ht

}

}

}

for (int i =0; i<2; i++)

{

for (int j =0; j<2; j++)

{

H[i][j] += A[i][j] - B[i][j]; //Теперьсчитаемновое Hi+1 = Hi+Ai-Bi

}

}

gradf0[0] = dfdx(x,0); //Подсчитаем градинет в т. Xi

gradf0[1] = dfdx(x,1);

p[0] = - (H[0][0]*gradf0[0]+H[0][1]*gradf0[1]);

p[1] = - (H[1][0]*gradf0[0]+H[1][1]*gradf0[1]); //Нашлиновое p = -H*f'(x)

}

cout<<" Х = [ "<<x[0]<<", "<<x[1]<<" ] "<<'\n';

cout<<"Значение аргумента F(X) = ";

cout<<fn(x);

getch();

return 0;

}

Поиск по сайту: