Newton-Raphson Method with Jacobian
27 次查看(过去 30 天)
显示 更早的评论
I have a problem with this program, a finite value vector is not returned despite the system having a solution.
Using function fsolve the result is Xeq3 = [0.6875 0.6346 0.9411], while using the function my_newton2 Xeq3 = [NaN NaN NaN].
I think the problem is in the function declaration but I was unable to find a solution. Please help me to find the error in the code.
function es()
KI = 1;
k1 = 1.2;
k2 = 1.1;
k3 = 1.3;
tf = 10;
X0 = [2, 1.25, 0.75];
%==================================
format long
syms x y z;
f(x,y,z) = [(k3 / (z + y) - k1 * x), (k1 * x - k2 * 0.75), (k2 * 0.75 - k3 * y)];
j(x,y,z) = jacobian(f, [x, y, z]);
j1(x,y,z) = inv(j);
Xeq3 = my_newton2(f, j1, [X0(1), X0(3), KI])
end
%===================================
function xval = my_newton2(f, j, x0)
tol = 1e-9; %tollerance
while true
% Calculate the next xval value
h = -f(x0(1),x0(2),x0(3))*j(x0(1),x0(2),x0(3));
xval = double(x0) + double(h);
%Check the result
exit = true;
for k=1:length(x0)
if abs(x0(k)-xval(k))>tol
exit=false;
end
end
if exit
break;
end
x0 = xval;
end
end
0 个评论
采纳的回答
J. Alex Lee
2020-7-21
In Newton loops you must evaluate your f and j and currently guessed iterate, so your line
h = -f(x0(1),x0(2),x0(3))*j(x0(1),x0(2),x0(3));
should rather look like (using your notation)
h = -f(xval(1),xval(2),xval(3))*j(xval(1),xval(2),xval(3));
so you will need a statement at the top to initialize xval to x0
5 个评论
J. Alex Lee
2020-7-21
If your newton's method works fine when your initial guess is close to the actual solution, then your newton's method is taking steps out of the region of convergence, and depending on your equations may be stepping into domains where your function values are NaN. This is just a limitation of newton's method, which requires reasonably good initial guesses. fsolve may be more robust to this issue.
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!