Newton-Raphson Method with Jacobian
32 views (last 30 days)
Show older comments
Morena Fabozzi
on 21 Jul 2020
Commented: Morena Fabozzi
on 21 Jul 2020
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 Comments
Accepted Answer
J. Alex Lee
on 21 Jul 2020
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 Comments
J. Alex Lee
on 21 Jul 2020
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.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!