Newtons method for system of nonlinear equations

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.

[Walter Roberson]

When the code invokes those cell arrays, what are you expecting will be returned?

{@(x,y,z)x+y+z-3;x^2+y^2+z^2-5;exp(x)+x*y-x*z-1}

means that you want to construct a cell array. The first element of the cell array is to be the handle to an anonymous function that expects three inputs and evaluates the sum of the inputs minus 3. The second element of the cell array looks in the current workspace for variables x, y, and z, in order to evaluate x^2+y^2+z^2-5 and drop the numeric result of the calculation into the cell array; this would fail if any of x, y, or z did not exist in the workspace that the cell array was being constructed in. The third element of the cell array looks in the current workspace for variables x, y, and z, in order to evaluate exp(x)+x*y-x*z-1 and drop the numeric result of the calculation into the cell array; this would fail if any of x, y, or z did not exist in the workspace that the cell array was being constructed in.

If all parts succeeded then the first parameter, f, would receive a cell array in which the first element was an anonymous function expecting 3 inputs, and the other two elements of the cell array were numeric values.

Then in the function, you try to index the cell array with x0, which contains [1,0,1] but 0 is an invalid index. If the x0 happened to contain only values that were 1, 2, or 3, then the indexing would work, returning a cell array, which you would then try to operate on with the * operator, which is not defined for cell arrays.

[Peter M Hogan]

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]

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.

Answer 1

[Walter Roberson]

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 .

[Peter M Hogan]

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

Answer 2

[Muhammad Asif]

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)

Answer 3

[ahmad haggag]

what is the final working code ?

[Peter M Hogan]

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