3.2.4. Розв'язування систем нелінійних рівнянь в matlab
Для розв'язання систем нелінійних рівнянь вигляду F(x) = 0 в MATLAB призначена функція fsolve. Тут x – вектор і F(x) – функція, яка повертає значення вектора.
Детальна інформація про функцію fsolve наведена в довідковій системі MATLAB. Наведемо лише короткий опис.
Синтаксис:
x=fsolve(fun, x0)
x=fsolve(fun, x0, options)
[x, fval, exitflag]=fsolve(...)
[x, fval, exitflag, output]=fsolve(...)
[x,fval,exitflag,output jacobian]=fsolve(...)
Опис:
fsolve знаходить корені (нулі) системи нелінійних рівнянь.
x = fsolve(fun, x0) – починає з точки x0 і намагається розв'язати описані в fun рівняння.
x = fsolve(fun, x0, options) – проводить оптимізацію з певними параметрами у структурній опції.
[x, fval, exitflag] = fsolve(...)– повертає значення exitflag, яке містить опис вихідних умов.
[x, fval, exitflag, output] = fsolve(...) – повертає структурний вихід з інформацією про оптимізацію
[x,fval,exitflag,output,jacobian] = fsolve(...)– повертає Якобіана від fun як розв'язок від х.
Приклад виклику функції:
x = fsolve(@myfun,x0) % myfun – функція користувача, що описана в m-файлі
x = fsolve(inline('sin(x.*x)'),x0);
Приклад. Нехай необхідно знайти розв'язок системи двох нелінійних рівнянь:
Створимо M-файл для розрахунку F(x).
function F = myfun(x)
F = [2*x(1) - x(2) - exp(-x(1));
-x(1) + 2*x(2) - exp(-x(2))];
За початкове наближення візьмемо точку x0=[-5; -5].
Викличемо підпрограму оптимізації
x0 = [-5; -5]; % початкове наближення
options=optimset('Display','iter'); %Опція відображення процесу обчислень
[x,fval] = fsolve(@myfun,x0,options) % виклик оптимізатора
Після виконання 33 ітерацій тримаємо розв'язок
x =
0.5671
0.5671
fval =
1.0e-006 *
-0.4059
-0.4059
Приклад знаходження розв'язку системи двох нелінійних рівнянь методами простих ітерацій та Ньютона в Matlab.
Нехай необхідно знайти розв'язок системи нелінійних рівнянь вигляду:
Створимо M-файл funF.m із заданим вхідним вектором функцій F:
function[F]=funF(x)
F=[f1(x(1),x(2)); f2(x(1),x(2))]; %тут необхідно задати вирази функцій f1 i f2
Наприклад, для системи
файл funF.m матиме вигляд:
function[F]=funF(x)
F=[x(1)^2+x(2)^2-4; x(1)*x(2)-1];
Також створимо окремий файл funJ.m для задання матриці Якобі:
function[J]=funJ(x)
J=[
;
% тут
потрібно задати елементи матриці Якобі,
які обчислюють аналітично
].
У нашому випадку файл funJ.m
function[J]=funJ(x)
J=[2*x(1), 2*x(2); x(2), x(1)] ; % тут потрібно задати елементи матриці Якобі
Тоді створимо M-файл mpi.m для реалізації методу простих ітерацій.
function [xout, dxout, mout]=mpi(x0, eps)
F=funF(x0); nf=norm(F); % знаходимо значення і норму вектора-функції в точці x0
J=funJ(x0); % знаходимо значення якобіана в точці x0
alpha=J(2,2)/(-J(1,1)*J(2,2)+J(1,2)*J(2,1));
beta=-alpha*J(1,2)/J(2,2);
gamma=-J(2,1)/(J(1,2)*J(2,1)-J(1,1)*J(2,2));
delta=-gamma*J(1,1)/J(2,1);
dx=[alpha*F(1)+beta*F(2); gamma*F(1)+delta*F(2)];
ndx=norm(dx);
m=0; x=x0;
xout=x'; dxout=dx'; mout=m;
while max([nf; ndx])>eps % перевіряємо виконання критерію закінчення ітерацій
x=x+dx; % знаходимо нове значення вектора розв'язку системи
F=funF(x); nf=norm(F); знаходимо значення і норму
dx=[alpha*F(1)+beta*F(2); gamma*F(1)+delta*F(2)]; % знаходимо чергову похибку
ndx=norm(dx); % знаходимо норму чергової похибки
m=m+1; % збільшуємо кількість ітерацій
xout=[xout; x']; dxout=[dxout; dx']; mout=[mout;m];
end;
Створимо M-файл kmn.m для реалізації класичного методу Ньютона.
function[xout, dxout, mout]=kmn(x0, eps)
F=funF(x0); nf=norm(F); % знаходимо значення функції і норму вектора-функції в точці x0
J=funJ(x0); % знаходимо значення якобіана в точці x0
delta=[det([F(:) J(:,2)]); det([J(:,1) F(:)])];
dx=-delta/det(J); ndx=norm(dx);
m=0; x=x0;
xout=x'; dxout=dx'; mout=m;
while max([nf; ndx])>eps % перевіряємо виконання критерію закінчення ітерацій
x=x+dx; % знаходимо нове значення вектора розв'язку системи
J=funJ(x);
F=funF(x); nf=norm(F);
delta=[det([F(:) J(:,2)]); det([J(:,1) F(:)])];
dx=-delta/det(J); ndx=norm(dx);
m=m+1;
xout=[xout; x']; dxout=[dxout; dx']; mout=[mout;m];
end;
Тоді у робочому середовищі Matlab задамо початкове наближення та точність обчислень:
x0=[2; 0.1];
eps=0.001;
Знайдемо розв'язок заданої системи рівнянь методом простих ітерацій:
[xout, dxout, mout]=mpi(x0, eps);
plot(mout, xout(:,1), mout, xout(:,2))
title('x mpi')
pause
plot(mout, dxout(:,1), mout, xout(:,2))
title('dx mpi')
pause
out_mpi=[mout, xout, dxout]
Отримаємо:
out_mpi =
0 2.0000 0.1000 -0.0226 0.4011
1.0000 1.9774 0.5011 -0.0407 0.0066
2.0000 1.9368 0.5077 -0.0026 0.0085
3.0000 1.9341 0.5162 -0.0019 0.0009
4.0000 1.9323 0.5171 -0.0003 0.0004
5.0000 1.9320 0.5175 -0.0001 0.0001
Обчислимо нев'язки:
>> funF([1.932, 0.518])
ans =
1.0e-003 *
0.9480
0.7760
Як бачимо, розв'язок x=[1.932, 0.518] з точністю до 0,001 знайдено за 5 ітерацій.
Знайдемо наближений розв'язок заданої СНР класичним методом Ньютона:
[xout, dxout, mout]=kmn(x0, eps);
plot(mout, xout(:,1), mout, xout(:,2))
title('x kmn')
pause
plot(mout, dxout(:,1), mout, xout(:,2))
title('dx kmn')
pause
out_kmn=[mout, xout, dxout]
Отримаємо:
out_kmn =
0 2.0000 0.1000 -0.0226 0.4011
1.0000 1.9774 0.5011 -0.0449 0.0159
2.0000 1.9326 0.5171 -0.0007 0.0006
3.0000 1.9319 0.5176 -0.0000 0.0000
Як бачимо, наближений розв'язок за класичним методом Ньютона знайдено за 3 ітерації.