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

3 次查看(过去 30 天)
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
Jan 2022-1-6
You forgot to mention, what the problem is. Why do you assume that there is something wrong?
Matthew Hunt
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.

请先登录,再进行评论。

采纳的回答

Jan
Jan 2022-1-11
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 个评论
Matthew Hunt
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.

请先登录,再进行评论。

更多回答(0 个)

产品


版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by