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

Question

Matthew Hunt 2022-1-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1623960-my-implemetation-for-newton-s-method-doesn-t-seem-to-be-working

评论： Matthew Hunt 2022-1-11

采纳的回答： Jan

在 MATLAB Online 中打开

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 个评论
显示无隐藏无

Jan 2022-1-6

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

Matthew Hunt 2022-1-7

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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Jan 2022-1-11

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1623960-my-implemetation-for-newton-s-method-doesn-t-seem-to-be-working#answer_872690

在 MATLAB Online 中打开

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Matthew Hunt 2022-1-11

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.

请先登录，再进行评论。