Newtons method for system of nonlinear equations

96 次查看(过去 30 天)
function p = sysNewton(f,J,x0,tol)
% f is the system of equations as a column vector
% this an anonymous function with a vector input and vector output
% J is the Jacobian of the system
% this is an anonymous function with a vector input and matrix output
% x0 is a set of initial guesses (in a column vector)
% tol is a tolerance
xold=x0;
xnew=x0-J(x0)^(-1)*(f(x0));
while norm(xnew-xold)>tol
xold=xnew;
xnew=xold-J(xold)^(-1)*f(xold);
end
p=xnew;
end
%code to call function
sysNewton({@(x,y,z)x+y+z-3;x^2+y^2+z^2-5;exp(x)+x*y-x*z-1},({@(x,y,z)1,1,1;2*x,2*y,2*z;y-z+exp(x),x,-x}),[1,0,1],10^-6)
My code is correct however I cannot call the function the way that I want to. I am misunderstanding how to call the function, as I didnt really need to know how to call it to write the script.
  3 个评论
Peter M Hogan
Peter M Hogan 2019-3-1
thanks walter. here is how i am calling my function now,
F =
function_handle with value:
@(x)[x(1)*exp(x(2))+x(3)-10;x(1)*exp(2*x(2))+2*x(3)-12;x(1)*exp(3*x(2))+3*x(3)-15]
>> k
k =
function_handle with value:
@(x)[exp(x(2)),x(1)*exp(x(2)),1;exp(2*x(2)),2*x(1)*exp(2*x(2)),2;exp(3*x(2)),3*x(1)*exp(3*x(2)),3]
>> p
p =
25 25 25
>> sysNewton(F,k,p,10^-6)
When I input this i get a 3x3 array. I do not understand what my answer is. What are the roots of the function?
Peter M Hogan
Peter M Hogan 2019-3-1
also my function is INCREDIBLY sensitive to intital guesses. I dont know I am lost and I have been working on this for way too long.

请先登录,再进行评论。

回答(3 个)

Walter Roberson
Walter Roberson 2019-3-1
编辑:Walter Roberson 2019-3-1
Your J(x0) returns 3 x 3 because of the way you define J. Raising that to power -1 using the ^ operator is inv() but more numerically unstable, and you should always avoid using inv() . You algebraic-matrix-multiply the inverse you calculate by the 3 x 1 matrix, giving a 3 x 1 result. If that was what you wanted then you would be better off coding J(x0)^(-1)*(f(x0)) as J(x0)\f(x0)
So you have a 3 x 1 on the right hand side. And you subtract it from x0, which is 1 x 3. In R2016a and earlier, that would be an error, but in R2016b and later, the effect is as if you had coded
bsxfun(@minus, x0, J(x0)etc)
giving a 3 x 3 result.
Your function will be **less* sensitive to initial guesses if you use the \ operator... but only relatively speaking.
By the way: my calculation is that there is no solution to that system, unless perhaps involving infinity. The first part, x(1) - FirstElementOf(J(x0)\F(x)) gives an equation that is 0 when x(2) = -2*log(2) or when x(2) is infinite. When you substitute in -2*log(2) then the third part, x(3) - ThirdElementOf(J(x0)\F(x)) gives 14/3 .
  1 个评论
Peter M Hogan
Peter M Hogan 2019-3-1
ok so what should i do instead? my jacobian should be 3x3 but I want my ouput "p" to be the roots of the system of the input function. which would be 3x1

请先登录,再进行评论。


Muhammad Asif
Muhammad Asif 2020-8-17
sysNewton(@(x) [x(1)+x(2)+x(3)-3;x(1)^2+x(2)^2+x(3)^2-5;exp(x(1))+x(1)*x(2)-x(1)*x(3)-1],...
@(x) [1,1,1;2*x(1),2*x(2),2*x(3);x(2)-x(3)+exp(x(1)),x(1),-x(1)],[1;0;1],10^-6)

ahmad haggag
ahmad haggag 2021-1-17
what is the final working code ?
  1 个评论
Peter M Hogan
Peter M Hogan 2021-1-17
function p = sysNewton(f,J,x0,tol)
% f is the system of equations as a column vector
% this an anonymous function with a vector input and vector output
% J is the Jacobian of the system
% this is an anonymous function with a vector input and matrix output
% x0 is a set of initial guesses (in a column vector)
% tol is a tolerance
xold=x0;
xnew=x0-J(x0)\(f(x0));
while norm(xnew-xold)>tol
xold=xnew;
xnew=xold-J(xold)\f(xold);
end
p=xnew;
end

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Operators and Elementary Operations 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by