solving numerically non-linear equation , code below does not work

algorithm to solve numrically set of two non-linear equations

%define functions , calculate Jacobian matrix ,  initial guess 
syms x y 
f(x,y) = x^2 + y - 2  ; % first function
g(x,y) = x + y^2 - 2  ; % second function
F_mat = [f ; g] ; % the RHS matrix in the algorithm
J(x,y)=jacobian([x^2 + y - 2 , x + y^2 - 2],[x,y]) ; % calculation of the jacobian matrix
%initialize the iterative process 
X0 = [1 ; 0] ; % initial guess vector
Xn = X0 ; % Xn is the vector x and y values that will zero the functions f and g
%for the iterative method 
    max_iter = 100 ; % max number of iterations
    tol = [1e-6 ; 1e-6] ; % tolerance for achieving the numerical solution 
    
    %loop of iterations till we converge to solution 
    for i = 1:max_iter
        F_mat = F_mat(X0) % i dont know how to calculate my matrix F_mat at my initial guess X0
        J = J(X0) % i dont know how to calculate my Jacobian matrix J at my initial guess X0
        Xn = X0 - J(X0)^-1 * F_mat(X0) % the Newton algorithm
        if abs(Xn-X0) < tol 
            break
        end
        
        Xn=X0 ;
    end

Error using indexing (line 249)
Invalid argument at position 2. Symbolic function expected 2 input arguments but received 1.

    %display X0
    disp(X0)

Answer 1

[John D'Errico]

What does not work? This is just a basic Newton-Raphson method. First, what are the functions for your test problem?

x^2 + y - 2

x + y^2 - 2

Before I do ANYTHING, what is the solution? Does a solution exist? If it does, the solution will be where the red and blue lines cross in this plot.

syms x y
fimplicit(x^2 + y - 2,'r')
hold on
fimplicit(x + y^2 - 2,'b')
hold off
grid on

Clearly, 4 solutions exist.

xy = solve(x^2 + y - 2 == 0,x + y^2 - 2 == 0)

xy = struct with fields:
    x: [4x1 sym]
    y: [4x1 sym]

[xy.x,xy.y]

ans = 

double([xy.x,xy.y])

ans = 4×2

1.0000 1.0000 -2.0000 -2.0000 1.6180 -0.6180 -0.6180 1.6180

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

So there are 4 solutions. A newton method will find one of them, IF it converges, and IF you did not start the solver out at a location where the solver diverges or gets stuck! My guess would be, if I start a solver out at [1,0], I'll end up at the first, or the third solution listed. But I might be wrong.

So next, write a basic Newton-Raphson scheme.

J = jacobian([x^2 + y - 2,x + y^2 - 2],[x,y])

J = 

Jfun = matlabFunction(J)

Jfun = function_handle with value:
    @(x,y)reshape([x.*2.0,1.0,1.0,y.*2.0],[2,2])

I've turned it into a function handle, so I can use it in a numerical method.

XY0 = [1;0];
obj = @(x,y) [x.^2 + y - 2;x + y.^2 - 2];

Note that both of these functions are a function of x and y, but I want to work with vectors. You can see how I did it below.

tol = 1e-8;
XY = XY0;
iter = 0;
itmax = 20;
deltaXY = inf;
while (deltaXY > tol) && (iter < itmax)
    iter = iter + 1;
    XYnew = [XY - Jfun(XY(1),XY(2))\obj(XY(1),XY(2))];
    deltaXY = norm(XYnew - XY);
    
    disp("Iteration: " + iter + ', DeltaXY: ' + deltaXY)
    
    XY = XYnew;
end

Iteration: 1, DeltaXY: 1.4142
Iteration: 2, DeltaXY: 0.4714
Iteration: 3, DeltaXY: 0.067344
Iteration: 4, DeltaXY: 0.0014328
Iteration: 5, DeltaXY: 6.4923e-07
Iteration: 6, DeltaXY: 1.3338e-13

format long g
XY

XY = 2×1

1.61803398874989 -0.618033988749895

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

It has converged quite rapidly to the third solution. Note my use of backslash on the jacobian, which implicitly does the equivalent of inv(Jfun(x,y))*obj(x,y),

[Yossi Mendel]

thanks , very impressive !

Appreciate if you tell me whats wrong with my code .

for example if i would like to calculate F_mat at x0 how do i write that ? F_mat(x0) ?

[Torsten]

%define functions , calculate Jacobian matrix ,  initial guess 
syms x y 
f = x^2 + y - 2  ; % first function
g = x + y^2 - 2  ; % second function
F_mat = [f ; g] ; % the RHS matrix in the algorithm
J=jacobian(F_mat,[x,y]);  % calculation of the jacobian matrix
%initialize the iterative process 
X0 = [1 ; 0] ; % initial guess vector
%for the iterative method 
max_iter = 100 ; % max number of iterations
tol = [1e-6 ; 1e-6] ; % tolerance for achieving the numerical solution 
%loop of iterations till we converge to solution 
for i = 1:max_iter
    F_mat_num = double(subs(F_mat,[x;y],X0)); % i dont know how to calculate my matrix F_mat at my initial guess X0
    J_num = double(subs(J,[x;y],X0)); % i dont know how to calculate my Jacobian matrix J at my initial guess X0
    Xn = X0 - inv(J_num) * F_mat_num; % the Newton algorithm
    if abs(Xn-X0) < tol 
        break
    end    
    X0=Xn ;
end
%display Xn
format long
disp(Xn)

   1.618033988749989
  -0.618033988749989

[Torsten]

Or nearer to your code:

%define functions , calculate Jacobian matrix ,  initial guess 
syms x y 
f(x,y) = x^2 + y - 2  ; % first function
g(x,y) = x + y^2 - 2  ; % second function
F_mat(x,y) = [f(x,y) ; g(x,y)] ; % the RHS matrix in the algorithm
J(x,y)=jacobian([x^2 + y - 2 , x + y^2 - 2],[x,y]) ; % calculation of the jacobian matrix

%initialize the iterative process 
X0 = [1 ; 0] ; % initial guess vector

%for the iterative method 
    max_iter = 100 ; % max number of iterations
    tol = [1e-6 ; 1e-6] ; % tolerance for achieving the numerical solution 
    
    %loop of iterations till we converge to solution 
    for i = 1:max_iter
        F_mat_iter = F_mat(X0(1),X0(2)) % i dont know how to calculate my matrix F_mat at my initial guess X0
        J_iter = J(X0(1),X0(2)) % i dont know how to calculate my Jacobian matrix J at my initial guess X0
        Xn = X0 - J_iter^-1 * F_mat_iter % the Newton algorithm
        if abs(Xn-X0) < tol 
            break
        end
        
        X0=Xn ;
    end

F_mat_iter = 

J_iter = 

Xn = 

F_mat_iter = 

J_iter = 

Xn = 

●

F_mat_iter = 

J_iter = 

Xn = 

F_mat_iter = 

J_iter = 

Xn = 

F_mat_iter = 

J_iter = 

Xn = 

    %display Xn
    format long
    disp(double(Xn))

   1.618033988749989
  -0.618033988749989

[John D'Errico]

How would I calculate F_mat? There are many ways we could do it, I am sure.

syms x y 
f(x,y) = x^2 + y - 2  ; % first function
g(x,y) = x + y^2 - 2  ; % second function

so we have f and g, as functions of x and y.

f_fun = matlabFunction(f,[x,y])

f_fun = function_handle with value:
    @(x,y)deal(y+x.^2-2.0,[x,y])

g_fun = matlabFunction(g,[x,y])

g_fun = function_handle with value:
    @(x,y)deal(x+y.^2-2.0,[x,y])

So we have two functions.

F_mat = @(x,y) [f(x,y);g(x,y)];

F_mat(2,3)

ans = 

I recall that [1,1] was one solution.

F_mat(1,1)

ans = 

[Yossi Mendel]

thanks much !