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

Рис. IV.4. Поверхность и карта линий уровня функции f x, y x2 y 2

Рис. IV.5. Поверхность и карта линий уровня функции f x, y x2 1.5 y 2

78

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

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

z=zeros(N);

i=1:N;

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

for j=l:N

z(i, j)=2*x(i); end;

end;

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

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

Z=zeros(N);

i=1:N;

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

z(i,j)=2*mu*x(i);

end;

end;

5.Создайте файл L_Grad.m (листинг IV.4), содержащий описание функции, возвращающей длину градиента функции (x,y,μ).

Листинг IV.4. Файл L_Grad.m function z=L_Grad(x,y,mu)

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

for j=l:N

z(i,j)=Dx_F_L4(x(i),y(j),mu) .^2+Dу_F_L4(х{i),у(j),mu) .^2; z(i,j)=z(i,j) .^0.5;

end;

end;

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

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

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

for i=1:N for j=1:N

z(i, j)=-Dx_F_L4(x(i), y(j), mu) ./L_Grad(x(i), y(j), mu); end;

end;

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

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

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

80

for j=l:N

z(i, j)=-Dy_F_L4(x(i), y(j), mu) ./L_Grad(x(i), y(j), mu);

end;

end;

8.Создайте файл Lambda.m (листинг IV.7), содержащий описание функции, возвращающей значение шага.

Листинг IV.7. Файл Lambda.m

function z=Lambda(x,у, nu, alpha, beta, gamma, lambda0); z=alpha*lambda0/ (beta+gamna*nu) ;

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

Листинг IV.8. Файл MG.m

function [X,Y,ff]=MG(x0,y0,mu,nu,alpha, beta, gamma, lambda0) X(l)=x0;

Y(1)=y0; ff(l)=F_L4(x0, y0, mu); for i=2:nu

X(i)=X(i-l) + Lambda(X(i-l),Y(i-1), nu, alpha, beta, gamma, lambda0)*s_x(X(i-1),Y(i-1),mu);

Y(i)=Y(i-l) + Lambda(X(i-l),Y(i-1), nu, alpha, beta, gamma, lambda0)*s_y (X(i-1),Y(i-1), mu);

ff(i)=F_L4(X(i),Y(j), mu);

end;

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

>> Vmax=20; % максимальное число шагов итерационного

%процесса

>>х0=2; у0=-1; % начальное приближение

>>lambda0=0.3; % начальное значение шага к экстремуму

81

>>beta=l; alpha=l; garama=l; % параметры, используемые

%для определения шага

>>mu=1;nu=20;

>>[X,Y,ff]=MG(x0,y0,mu.nu,alpha,beta.gamma,lambda0);

>>plot(X);

>>hold on

>>plot(Y); % номера итерационного процесса (рис. IV.6)

>>hold off;

>>plot(ff); % зависимость значения функции f(x,y) от номера

%итерации (рис. IV.7)

>> stairs(X,Y); % построение траектории итерационного

%процесс;

%на плоскости ХоY (рис. IV.8) .

>>plot3(X,Y,ff) % построение траектории итерационного

%процесса Б трехмерном пространстве(рис. IV.9)

Рис. IV.6. Траектория решения системы нелинейных

уравнений на плоскости ХоY

82

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

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

Рис. IV.9. Визуализация итерационного процесса в трехмерном пространстве

Пример IV.2. Алгоритм поиска экстремума с шагом, зависящим от свойств минимизируемой функции (использование производной по направлению). Исследуем данный алгоритм применительно к минимизации функции двух переменных, заданной полиномом 4-го порядка:

f x, y kx x4 ky y4 lx x3 ly y3 ox x2 oy y2 px x py y q kxy x4 y4 lxy x3 y3 oxy x2 y2 pxy xy

Форма функции определяется коэффициентами kx, ky, lx, ly ит.д. 1. Создайте файл F2_L4.m (листинг IV.9), содержащий описание функции, возвращающей значения функции/(х, у).

Листинг 4.9. Файл F2_L4.m

functionz=F2_L4(х,у,Lx,Lу,Ох,Рх,Lxy,Оху,O,Кх,Ку,Oу,Ру,Кху,Рху) N=length(x);

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

for j=l:N

z(i,j)=Kx*x(i) .^4+Ky*y(j) .^4+Lx*x(i) .^3+Ly*y(j) .^3+...

Ox*x(i) .^2+Qy*y(j) .^2+Px*x(i)+Py*y(j)+Q+...

84

Kxy*x(i) .^4*y(j) .^4+Lxy*x(i) .^3*y(j) .^3+...

Oxy*x(i) .^2*y(j) .^2+Pxy*x(i) .*y(j) ;

end;

end;

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

»Lx=0.3;Ly=0.4;Ох=1;Рх=1.2;Lxy=0.4;Оху=0.3;Q=3;

»Кх=0.5;Ку=1;0у=2;Ру=0.6;Кху=0.2;Рху=0.1;

»Xmin=-l;Xmax=l;

»Ymin=-l;Ymax=l;

»N=23;

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

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

»M=F2_L4(x,y,Lx,Ly,Ox,Px,Lxy,Oxy,Q,Kx,Ky,Oy,Py,Kxy,Pxy) ;

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

»surfc(X,Y,z)

Рис. IV.10. Визуализация функции f (x, у)

— График функции, изображенной на рис. IV.10, имеет "плато" — чрезвычайно пологое "дно". Понятно, что градиент в

85

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

3. Создайте файлы g2_x.m (листинг IV.10) и g2_y.m (листинг IV.11), содержащие описание функций, возвращающих значения частных производных.

Листинг IV.10. Файлы g2_x.m

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

z=zero9(N); for i=l:N for j=l:N

z(i, j)=4*Kx*x(i) .^3+3*Lx*x{i) .^2+2*Ох*х(i)+Px+.. . 4*Kxy*x(i) .^3*y(j).^4+3*Lxy*x(i) .^2*y(j) .^3+...

2*0xy*x(i) .*y(j) .^2+Pxy*y(j);

end;

end;

Листинг IV.11. Файлы g2_y.m function

z=g2_y(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)=4*Ky*y(j).^3+3*Ly*y(j).^2+2*Oy*y(j)+...

Py+4*Kxy*x(i) .^4*y(j) .^3+3*Lxy*x(i) .^3*y(j) .^2+...

2*Oxy*x(i).^2*y(j)+Pxy*x(i);.

end;

end;

86

4.Создайте файл L2_Grad,m (листинг IV.12), содержащий описание скалярной функции, возвращающей значение длины градиента функции f(x,y).

Листинг IV.12. Файл L2_Grad.m

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

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

for j=l:

z(i,j)= L2_Grad (x(i) ,y(j) ,Lx,Ly,Ox, Px,Lxy,Oxy,Q,Kx,Ky,Oy,Py,...

Kxy.Pxy).^2+...

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

z(i,j)=z(i,j) .^0.5; end;

end;

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

Листинг IV.13. Файлы S2_x.m

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

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

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

for j=1:N

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

end;

end;