Материал: Спецглавы высшей математики. численные методы. Пантелеев И.Н

Листинг IV.14. Файлы S2_y.m

function z=S2_y(x,y,Lx,Ly,0x,Px,Lxy,Oxy,Q,Kx,Ky,Oy,Py,Kxy,Pxy); N=length(x) ;

z=zeros(N);

i=1:N;

j=1:N; for i=l:N

for j=l:N

z(i,j)=-g2_y(x(i),y(j),Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Oy,Py,Kxy,Pxy)./...

. L2_Grad(x(i),y(j),Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Oy,Py,Kxy,Pxy); end;

end;

6. Создайте файлы g2_xx.m. (листинг IV.15), g2_xy.m (листинг IV.16), g2_yy.m (листинг IV.17), содержащие описание функ-

Листинг IV.15. Файлы g2_xx.m

function z=g2_xx(x,y,Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Oy,Py,Kxy,Pxy) N=length(x);

z=zeros(N); for i=l:N

for j=l:N z(I,j)=12*Kx*x(i).^2+6*Lx*x(i)+2*Ox…

+12*Kxy*x(i).^2*y(j).^4+6*Lxy*x(i).*y(j).^3+2*Oxy*y(j).^2; end;

end;

Листинг IV.16. Файлы g2_xy.m

function z=g2_xy(x,y,Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Oy,Py,Kxy,Pxy) N=length(x);

z=zeros(N); for i=1:N

for j=l:N

88

z(i,j)=16*Kxy*x(i).^3*y(j).^3+9*Lxy*x(i).^2*y(j).^2+… 4*Oxy*x(i) .*y(j)+Pxy;

end;

end;

Листинг IV.17. Файлы g2_yy.m function

z=g2_уу(х,у,Lх,Lу,Ох,Рх,Lху,Оху,С,Кх,Ку,Оу,Ру,Кху,Рху) N=length(x);

z=zeros(M); for i=1:N for j=l:N

z(i,j)=12*Ky*y(j).^2+6*Ly*y(j)+2*Oy+...

12*Kxy*x(i).^4*у(j).^2+6*Lху*х(i).^3*у(j)+2*Oxy*x(i).^2; end;

end;

7. Создайте файл Lambda2.m (листинг IV.18), содержащий описание функции, возвращающей значения коэффициента λ.

Листинг IV.18. Файл Lambda2.m

function z=Lambda2(x,y,Lx,Ly,Ox,Px,Lxy,Oxy,Q.,Kx,Ky,Qy,Py,Kxy,Pxy) N=length(x);

z=zeros(N); for i=l:N

for j=l:N al=g2_x(x(i),y

(j),Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Oy,Py,Kxy,Pxy).*...

S2_x(x(i),y(j),Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Oy,Py,Kxy,Pxy); a2=g2_y(x(i),y(j),Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Oy,Py,Kxy,Pxy).*...

S2_y(x(i),y(j),Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Oy,Py,Kxy,Pxy); b1=g2_xx(x(i),y(j),Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Cy,Fy,Kxy,Pxy).* S2_x(x(i),y(j),Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Oy,Py,Kxy,Pxy).^2; b2=2*g2_ху(х(i),y(j),Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Qy,Py,Kxy,Pxy)*_ S2_x(x(i),y(j),Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Oy,Py,Kxy,Pxy).*...

S2_y(x(i),y(j),Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Oy,Py,BOty.Pty);

89

g3=g2_yy(x(i),y(j),Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Oy,Py,Kxy,Pxy).*...

S2_y(x(i),y(j),Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Oy,Py,Kxy,Pxy).^2; z(i.j)=-(a1+a2)./(b1+b2+b3);

end;

end;

8. Создайте файл MG2.m (листинг IV.19), содержащий описание функции, возвращающей значения переменных х, у и соответствующее значение функции f(x,y,μ) на каждом шаге итерационного процесса.

Листинг IV.19. Файл MG2.m

function [X,Y,ff]=MG(xO,yO,nu,Lx,Ly,Ox,Px,Lxy,Oxy,… ….Q,Kx,Ky,Cy,Py,Kxy,Pxy)

Х(1)=х0; Y(l)=y0;

ff(l)=F2_L4(x0,y0,Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Oy,Py,Kxy,Pxy); fori=2:nu

X(i)=X(i-l)+Lambda2(X(i-l),Y{i-l),Lx,Ly,Ox,Px,Lxy,Oxy,Q,._ Kx,Ky,Oy,Py,Kxy,Pxy)*s2_x(X(i-l),Y(i-l),Lх,Lу,Ох,Рх,Lху,Оху,

Q,Kx,Ky,Oy,Py,Kxy,Pxy), Y(i)=Y(i-l)+Lambda2(X(i-l),Y(i-l),Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,...

Ky,Oy,Py,Kxy,Pxy)*s2_y(X(i-l),Y(i-l),Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,...

Oy,Py,Kxy,Pxy); ff(i)=F2_L4(X(i),Y(i),Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Oy,Py,Kxy,Pxy);

end;

9. Найдите минимум функции f(x,y,μ) и визуализируйте итерационный процесс.

»хО=1;уО=0.5; % начальное приложение

»nu«10; % число шагов итерационного процесса

»[ХУff]=MG2(xO,yO,nu,Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Oy,PY, Kxy.Pxy); » plot(X); hold on; plot(Y); hold off % рис. IV.11

»plot(ff); % рис. IV.12

»stairs(X,Y); % рис. IV.13

»plot3(X,Y,ff); % рис. IV.14

Обратите внимание на чрезвычайно быструю сходимость и на отсутствие колебаний решения. Это достигается за счет того, что формула для шага строится с учетом формы минимизируемой функции.

Из графика, представленного на рис. IV.10, видно, что функция имеет очень пологое "дно", поэтому для алгоритмов с шагом, не зависящим от формы функции, это было бы причиной очень плохой сходимости.

Рис. IV.11. Зависимость координат точки решения от номера итерации

91

Рис. IV.12. Зависимость значения исследуемой функции от номера итерации

Рис. IV.13. Траектории итерационного процесса в плоскости XoY

92

Рис. IV.14. Траектория итерационного процесса

впространстве

4.4.Метод Ньютона

Метод Ньютона основан на квадратической аппроксимации минимизируемой функции в окрестности точки х(k). Минимум квадратической функции легко найти, приравнивая ее градиент к нулю. Можно сразу же вычислить положение экстремума и выбрать его в качестве следующего приближения к точке

F x k 1 F x k x k Fsq x k 1 ,

где

93

Для учета фактических особенностей минимизируемой

функции будем использовать в (IV.22) значения градиента и матрицы вторых производных, вычисленных не по аппроксимирующей Fsq(.), а непосредственно по минимизируемой функ-

ции F(x). Заменяя Fsq(.) в (4.22), найдем длину шага x k :

Метод Ньютона реализуется следующей последовательностью действий:

1.Задается (произвольно) точка начального приближения x 0 . 2.Вычисляется а цикле по номеру итерации k=0,1,:

94

• значение вектора градиента F x k в точкеx x k ;

•значение матрицы вторых производных 2 F x k ;

•значение матрицы обратной матрице вторых производных

2 F x k 1 ;

• значение шага x k по формуле (IV.23);

• новое значение приближения x k 1 по формуле (IV.19). 3.Закончить итерационный процесс, используя одно из условий, описанных в разд. 2.5.

Оценка погрешности метода Ньютона дается формулой

аналогичной формуле, полученной при доказательстве теоремы о достаточном условии сходимости итерационного процесса.

Достоинства метода Ньютона.

•Для квадратической функции метод позволяет найти минимум за один шаг.

•Для функций, относящихся к классу поверхностей вращения (т. е. обладающих симметрией), метод также обеспечивает сходимость за один шаг (поскольку в точке минимума аргументы минимизируемой функции и ее квадратической аппроксимации совпадают).

•Для несимметричных функций метод обеспечивает более высокую скорость сходимости, чем при использовании других модификаций метода наискорейшего спуска.

Кнедостаткам метода Ньютона следует отнести необходимость вычислений и (главное!) обращения матриц вторых производных. При этом не только расходуется машинное время,

95

но, что более важно, могут появиться значительные вычислительные погрешности, если матрица 2 F x k окажется плохо

обусловленной.

Примечание.

В качестве количественной характеристики обусловленности матрицы А, обозначаемой condA принимают произведение

A A 1 , где — норма матрицы, определяемая по аналогии с нормой вектора:

A 2 max n aik .

4 1 m

k 1

Можно показать, что при использовании нормы A1 чис-

ло обусловленности матрицы связано с максимальным собственным значением матрицы АА* (здесь А* — эрмитово сопря-

f x, y 100 y x2 2 1 x 2 методом, в котором шаг зависит

от свойств минимизируемой функции (метод Ньютона). Создайте файл Rozenbrok.m (листинг IV.20), содержащий описание функции, возвращающей значения функции Розенброка.

Листинг IV.20. Файл Rozenbrok.m function z=Rozenbrok(x,y) N=length(x) ;

z=zeros (N) ; for i=l:N

for j=l:N z(i,j)=100*(y(j)-x(i).^2).^2+(l-x(i)).^2; end;

end;

96

2. Задайте координатную сетку.

»N=23; %задание числа узлов координат

»i=1 :N; j=l:N;

»Xmin=-0.2; Xmax=1.2;

»Ymin=-0.2; Ymax=1.2;

% задание координат узлов сетки

»х(i)=Xmin+(Xmax-Xmin)/(N-l)*(i-1);

»y(j)=Ymin+(Ymax-Ymin)/(N-1)*(j-1);

3. Вычислите значения функции в узлах координатной сетки. » z=Rozenbrok(x,y) ;

4.Визуализируйте функцию Розенброка и карту линий уровня.

»[X Y]=meshgrid(x,y) ;

»surf(X,Y,Z); % визуализация поверхности, задаваемой функцией % Розенброка (рис. IV.15)

»contour (X,Y, z, 53) ; % визуализация карты эквипотенциалей

%функции Розенброка (рис. IV.16)

5.Создайте файл R_x.m (листинг IV.21), содержащий описание функции, возвращающей значение первой производной по переменной х.

Листинг IV.21. Файл R_x.m function z=R__x(x,y) N=length(x);

z=zeros(H); for i=sl;N

for j=l:N

z(i, j)=-400*(y(j)-x(i) .^2) .*x(i)-2*(l-x(i)); end;

end;

6.Создайте файл R_y.m (листинг IV.22), содержащий описание функции, возвращающей значение первой производной по переменной у.

97