Материал: Спецглавы высшей математики. численные методы. Пантелеев И.Н
Рие. IV.15. Поверхность, задаваемая функцией Розенброка
Рис. IV.16. Карта линий уровня функции Розенброка
98
Листинг IV.22. Файл R_y.m function z=R_y(x,y) N=length(x);
Z=zeros(N); for i=l:N
for j=l:N z(I,j)=200*(y(j)-x(i).^2);
end;
end;
7.Создайте файл R_xx.m (листинг IV.23), содержащий описание функции, возвращающей значение второй производной по переменной х.
Листинг IV.23. Файл R_xx.m function z=R_xx(x,y) N=length(x);
Z=zeros(N); for i=l:N for j=l:N
z(i,j)=1200*x(i).^2-400*y(j)+2; end;
end;
8.Создайте файл R_yy.m (листинг IV.24), содержащий описание функции, возвращающей значение второй производной по переменной у.
Листинг IV.24. Файл R_yy.m function z=R_yy(x,y) N=length(x);
z=zeros (N); for i=l:N
for j=l:N z(i,j)=200;
end;
end;
99
9. Создайте файл R_xy.m (листинг IV.25), содержащий описа- |
B(2,l)=R_y(x,y); |
|
ние функции, возвращающей значение смешанной производной |
z=A^-l*B; |
|
по переменным х, у. |
|
12.Визуализируйте итерационный процесс (рис. IV.17-IV.19). |
Листинг IV.25. Файл R_xy.m |
|
» plot(x,'-x'); |
function z=R_xy(x,y) |
|
» hold on |
N=length(x); |
|
» plot(у,'-о'); |
z=zeros(N); |
|
» hold off |
for i=l:N |
|
» plot(z,'-o') |
for j=l:N |
|
» plot3(x,y,z,'-o'); grid on |
z(i,j)=-400*x(i); |
|
|
end; |
|
|
end; |
|
|
10. Создайте файл Rg.m (листинг IV.26), содержащий описание |
|
|
функции, возвращающей значения обратной матрицы вторых |
|
|
производных на вектор-градиент. |
|
|
Листинг IV.26. Файл Rg.m |
|
|
function z=Rg(x,y) |
|
|
A(l,l)=R_xx(x,y); |
|
|
A(l,2)=R_xy(x,y); |
|
|
A(2,l)=R_xy(x,y); |
|
|
A(2,2)=R_yy(x,y); |
|
|
B(l,l)=R_x(x,y); |
|
|
B(2,l)=R_y(x,y); |
|
|
z=А^-1*В; |
|
|
11.Создайте файл RnO.m (листинг IV.27), содержащий описа- |
Рис. IV.17. Зависимость значений координат решения уравне- |
|
ние функции, возвращающей значения переменных и соответ- |
ния от номера итерации |
|
ствующих значений функции Розенброка. |
|
|
Листинг IV.27. Файл RnO.m |
|
|
function z=Rg(x,y) |
|
|
A(l,l)=R_xx(x,y); |
|
|
A(l,2)=R_xy(x,y); |
|
|
A(2,l)=R_xy(x,y); |
|
|
A(2,2)=R_yy(x,y); |
|
|
B(l,l)=R_x(x,y); |
100 |
101 |
|
||
Рис. IV.18. Зависимость значения исследуемой функции от номера итерации
Рис. IV.19. Траектория итерационного процесса в пространстве
102
4.5. Решение систем нелинейных уравнений средствами пакета MATLAB
Средствами пакета MATLAB найдем решение системы нелинейных уравнений
x2 y 2 4 0
y x2 1 0
которая может быть записана в векторном виде
|
|
|
|
x 0, |
|
|
|||
F |
|
|
|||||||
где |
|
|
|
|
|
|
|||
|
|
|
|
|
2 |
2 |
|
|
|
F x |
x1 |
x2 |
4 |
|
|||||
x |
2 |
x2 |
1 |
, |
|||||
|
|
|
|
|
1 |
|
|
||
xx1 .
x2
Для этого необходимо выполнить следующую последовательность действий:
1. Создать файл fm.m (листинг IV.28), содержащий описание функции, возвращающей значения функции F x .
Листинг IV.28. Файл fm.m function z=fm(x) z(l,l)=x(l).^2+х(2)^2-4; z(2,l)=x(2)-x(l)^2-l;
2. Задать вектор начального приближения.
»z(1,1)=1;
»z(2,1)=1;
3.Обратиться к встроенной функции fsolve(), возвращающей решение системы линейных уравнений.
» х = fsolve ( 'fm' ,z,optimset( 'f solve' ) ) Optimization terminated successfully:
103
Relative function value changing by less than OPTIONS . TolFun
%Процесс оптимизации завершен
%Относительная величина изменения функции меньше чем
%переменная OPTIONS . TolFun
х = 0.8895 1.7913
4.Проверить полученное решение.
»fm(x) ans =
l.0e-008 * 0.5162 -0.4888
Для одновременного вывода координат вектора решения уравнения и соответствующего значения вектор-функции следует выполнить команду:
»[х fval] = f solve ( 'fm' ,z,optimset('f solve'))
Optimization terminated successfully:
Relative function value changing by less than OPTIONS . TolFun x =
0.8895
1.7913
fval = l.0e-008 * 0.5162 -0.4888
Для вывода на экран монитора значений вектора-решения и соответствующего значения функции F x следует выполнить
команду:
» [х fval exitflag] = fsolve('fm',z,optimset('Display','Iter'));
104
|
|
Norm of |
|
First-order |
|
|
Iteration |
Func-count f(x) step |
optimality |
CG-iterations |
|||
1 |
4 |
5 |
1 |
|
5 |
0 |
2 |
7 |
1 |
1 |
|
4 |
1 |
3 |
10 |
0.0026 |
0.223607 |
0.17 |
1 |
|
4 |
13 |
4.53078e-008 |
0.013546 |
0.00055 |
1 |
|
5 |
16 |
5.05378e-017 |
7.18274e-005 1.79e-008 |
1 |
||
Optimization terminated successfully:
Relative function value changing by less than OPTIONS.TolFun
В ряде случаев более удобным оказывается метод Ньютона, при использовании которого (как описано в разделе 4.2) необходимо знать в данной точке и значения функции, и значения якобиана, что позволяет реализовать итерационный процесс. Метод Ньютона в пакете MATLAB реализуется следующей последовательностью действий:
1. Создайте файл Fml.m (листинг IV.29), содержащий описание функции, возвращающей одновременно значения функции
F x и значения якобиана.
Листинг IV.29. Файл Fml.m function [z,J]=fm(x) z(l,l)=x(l).^2+х(2)^2-4; z(2,l)=x(2)-x(l)^2-l; J(l,l)=2*x(l);
J(l,2)=2*x(2);
J(2,l)=-2*x(l); J(2,2)=l;
2.Включите режим использования метода Ньютона и отображения итерационного процесса на экране.
» options=optimset('Jacobian','on','Display','Iter');
105
3.Обратитесь к встроенной функции fsolve(), возвращающей
решение системы линейных уравнений. |
|
|
|
|||||
» [х fval exitflag] = fsolve('fml',z,options); |
|
|
|
|||||
|
|
Normof First-order |
|
|
|
|
||
Iteration |
Func-count |
f(x) |
step |
optimality |
CG-iterations |
|||
1 |
2 |
|
5 |
1 |
|
5 |
|
0 |
2 |
3 |
|
1 |
1 |
|
4 |
|
1 |
3 |
4 |
|
0.0026 |
0.223607 |
|
0.17 |
|
1 |
4 |
5 |
4.53076e-008 |
0.013546 |
|
0.00055 |
1 |
||
5 |
6 |
5.04984e-017 |
7.18272e-005 |
1.79e-008 |
1 |
|||
Optimization terminated successfully:
Relative function value changing by less than OPTIONS.TolFun
Рассмотрим реализацию метода спуска средствами пакета MATLAB на примере системы нелинейных уравнений, который был разобран ранее в настоящем разделе. Напомним, что для этого (следуя подходу, описанному в разделе 4.3) нужно из уравнений исходной системы создать новую положительно определенную функцию, минимум которой и будет искомым решением системы нелинейных уравнений. Следовательно, для нахождения решения рассматриваемой системы нелинейных уравнений необходимо выполнить следующую последовательность действий:
1. Создать файл F_sq.m (листинг IV.29), содержащий описание функции, возвращающей значения суммы квадратов функций f(x, у) = х2 + у2 -4 и g(x,y)=y-x2-1.
Листинг IV.30. Файл F_sq.m function z=F_sq(x)
s=fm(x); z=s(1,1).^2+s(2,1).^2;
Здесь мы предполагаем, что ранее файл fm.m, содержащий описание функции, возвращающей значения функций f (x,y) и g(х, у), уже создан.
2.Задать начальное приближение, в окрестности которого будет находиться минимум функции f(х, у)2 + g(x, у)2 .
» x=[1;1];
3.Обратиться к встроенной функции fminsearch ()
» fminsearch ('f_sq',х) ans =
0.8896
1.7913
Способы обращения к функции fminsearch( ) и список ее формальных параметров аналогичны способам обращения к функции fsolve () и списку ее формальных параметров.
107
106