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

4 visualizzazioni (ultimi 30 giorni)
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 Commenti
Jan
Jan il 6 Gen 2022
You forgot to mention, what the problem is. Why do you assume that there is something wrong?
Matthew Hunt
Matthew Hunt il 7 Gen 2022
Oh, sorry, the answer diverges rather than converges. Even when you're close to the solution it tends to infinity.

Accedi per commentare.

Risposta accettata

Jan
Jan il 11 Gen 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 Commento
Matthew Hunt
Matthew Hunt il 11 Gen 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.

Accedi per commentare.

Più risposte (0)

Categorie

Scopri di più su Programming in Help Center e File Exchange

Prodotti


Release

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by