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

4 ビュー (過去 30 日間)
Matthew Hunt
Matthew Hunt 2022 年 1 月 6 日
コメント済み: Matthew Hunt 2022 年 1 月 11 日
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 件)

カテゴリ

Help Center および File ExchangeProgramming についてさらに検索

製品


リリース

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by