利用牛頓法接非線性方程組的Matlab程式例項
阿新 • • 發佈:2018-12-26
首先將線性方程組寫成f(x)=0的形式,編寫第一個Matlab函式
function f = fun()
% 求解:
% 3*x1-cos(x2*x3)-1/2 = 0
% x1^2-81*(x2+0.1)^2+sin(x3)+1.06 = 0
% exp(-x1*x2)+20*x3+(10*pi-3)/3 = 0
% 求解精度為0.00001
syms x1 x2 x3
f1 = 3*x1-cos(x2*x3)-1/2;
f2 = x1^2-81*(x2+0.1)^2+sin(x3)+1.06;
f3 = exp(-x1*x2)+20*x3+(10*pi-3)/3;
f = [f1 f2 f3];
然後是方程組的Jacobi矩陣,編寫第二個Matlab函式
function df = dfun()
f = fun();
df = [diff(f,'x1'); diff(f, 'x2'); diff(f, 'x3')];
df = df';
然後求解方程組的程式
function x = newton(x0, eps, N)
% x0是初值,可以任取,不需要滿足方程組
% 最後得到的解與初值選取有關
% eps為精度
% N最大迭代次數
for i = 1:N
f = eval(subs(fun(), {'x1', 'x2', 'x3'}, {x0(1), x0(2), x0(3)}));
df = eval(subs(dfun(), {'x1' , 'x2', 'x3'}, {x0(1), x0(2), x0(3)}));
x = x0 - f/df;
if norm(x - x0) < eps
for j = 1:length(x)
fprintf('%.2f\t', x(j));
end
break;
end
x0 = x;
end