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

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?

[Jan]

You forgot to mention, what the problem is. Why do you assume that there is something wrong?

[Matthew Hunt]

Oh, sorry, the answer diverges rather than converges. Even when you're close to the solution it tends to infinity.

Answer 1

[Jan]

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

[Matthew Hunt]

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.