My implemetation for Newton's method doesn't seem to be working

3 views (last 30 days)
Matthew Hunt
Matthew Hunt on 6 Jan 2022
Commented: Matthew Hunt on 11 Jan 2022
My goal is to solve the following set of nonlinear equations:
for this I use the Newton method in the form:
The code to define my equations is given by:
function Fn=f_test(x)
Fn=[x(1)^2-x(2)+x(3)^3-26;x(1)+x(2)+x(3)-6;x(1)^2+x(2)^3-x(3)^2];
The code to compute the Jacobian is given by:
function [J]=mat_jac(x,f)
% computes the Jacobian of a function
n=length(x);
eps=1e-5; % could be made better
J = zeros(n,n);
T=zeros(1,n);
for i=1:n
T(i)=1;
x_plus = x+eps*T;
x_minus = x-eps*T;
J(:,i) = (f_test(x_plus)-f_test(x_minus))/(2*eps);
T(i)=0;
end
My test script is given by:
%This tests the algorithm
x=[0.5 0.5 0.5];
x_old=x';
end_error=1e-5;
error=10;
while (error>end_error)
J=mat_jac(x_old',f_test(x_old'));
F=-f_test(x_old);
dx=linsolve(J,F);
x_new=x_old+dx;
x_old=x_new;
error=norm(f_test(x_old));
end
The solution to these equations is x=1,y=2,z=3 and I'm starting the search off at x=0.5,y=0.5,z=0.5. I know my jacobian is correct as I've checked it and my Newton's method looks correct so I don't know where I'm going wrong. Any suggestions?
  2 Comments
Matthew Hunt
Matthew Hunt on 7 Jan 2022
Oh, sorry, the answer diverges rather than converges. Even when you're close to the solution it tends to infinity.

Sign in to comment.

Accepted Answer

Jan
Jan on 11 Jan 2022
I've cleanup the code, e.g. removed the not needed "x_old". Avoid using the names of important Matlab functions as variables: eps, error. Avoid mixing row and column vectors frequently, but just work with column vectors in general.
Finally, it does converge.
x = [0.5; 0.5; 0.5];
limit = 1e-5;
err = 10;
loop = 0;
while err > limit
J = mat_jac(x, @f_test); % 2nd argument is not f_test(x)
F = f_test(x);
x = x - J \ F;
err = norm(f_test(x))
loop = loop + 1;
if loop == 1e5
error('Does not converge!');
end
end
err = 1.1332e+03
err = 334.4893
err = 135.1517
err = 97.7923
err = 52.6326
err = 317.6439
err = 7.0341e+04
err = 2.0854e+04
err = 6.1628e+03
err = 1.8123e+03
err = 526.0138
err = 147.4260
err = 37.6927
err = 7.6553
err = 0.9596
err = 0.0499
err = 2.4892e-04
err = 6.4461e-09
function J = mat_jac(x, f)
n = length(x);
v = 1e-8;
J = zeros(n, n);
T = zeros(n, 1);
for i=1:n
T(i) = 1;
J(:, i) = (f(x + v * T) - f(x - v * T)) / (2 * v);
T(i) = 0;
end
end
function Fn = f_test(x)
Fn = [x(1)^2 - x(2) + x(3)^3 - 26; ...
x(1) + x(2) + x(3) - 6; ...
x(1)^2 + x(2)^3 - x(3)^2];
end
  1 Comment
Matthew Hunt
Matthew Hunt on 11 Jan 2022
The strange thing that the next day when I tried it, it worked! I used as an initial guess for initially and it started working, and as I increased epsilon, it began to work and then when I went to the old settings again, it still worked. Very odd.
As for the x_old variable, I prefer it that way.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!