Newton's method for 2 dimension vectors

I want to write a code which gives me the roots of a vector function with sing the Newton's method.

that's what I wrote:

%%file newton2d.m
function [out] = newton2d (f,Jf,x,y, eps)
    k=0;
    Jf = jacobian(f, [x, y]) %error
    h = f(x,y)./Jf(x,y); %error
    while abs(h)< eps && (k < 1000)
        h = f(x,y)./Jf(x,y);
        [x,y] = [x,y] - h;
        k = k+1;
    end
    out = [x,y];
end

But there is a problem with the jacobian: i don't know where to insert it?

What can I do to make newton2d function work?

Thank you !

[Kevser Cifci]

I changed a little bit my code:

%%file newton2d.m
function [out] = newton2d(x, y, f, Jf, eps)
    k = 0;
    h = Jf(x,y)\f(x,y)';
    while abs(h)>=eps && k<1000
        h = Jf(x,y)\f(x,y)'; %% error
        [x,y] = [x,y]- h'; 
        k = k+1;
    end
    out = [x,y];
end

Why it's still not working?

Then I create a vector function with it's Jacobian:

%%file f4.m
function out = f4(x,y)
    out = [(x^3)-(3*x*(y^2))-1 , (3*(x^2)*y)-(y^3)]
end

%%file Jf4.m
function out = Jf4(x,y)
    out = jacobian(f, [x, y]); 
end

How to call f4 in my jacobian function?

newton2d(1,2,@f4, @Jf4, 10^(-8))

it's not working when I call my functions...

Thank you for your help !

Answer 1

[Alan Stevens]

LIke this?

% Functions
 f = @(XY) [XY(1).^3 - 3*XY(1).*XY(2).^2 - 1;
            3*XY(1).^2.*XY(2) - XY(2).^3];
 J = @(XY) [3*XY(1).^2 - 3*XY(2).^2,  -6*XY(1)*XY(2); 
            6*XY(1).*XY(2), 3*XY(1).^2 - 3*XY(2).^2];
% Initial guesses
XY = [2; 1];
tol = 10^-6;
err = 1;
maxits = 100;
its = 0;
while (err>tol) && (its<maxits)
     XYold = XY;
     XY = XYold - J(XYold)\f(XYold);
     err = norm(XY-XYold);
     its = its+1;
end
x = XY(1); y = XY(2);
disp('parameter values:')

parameter values:

disp(['x = ',num2str(x), '  y = ', num2str(y)])

x = 1  y = -1.0223e-17

disp('function values:' )

function values:

disp(f(XY))

   1.0e-16 *

         0
   -0.3067

[Kevser Cifci]

Yes that's the idea but it has to be a general newton2d function which can be applied to others vectors functions.

[Alan Stevens]

More like this then:

% Functions 
 f = @(xy) [xy(1).^3 - 3*xy(1).*xy(2).^2 - 1;
            3*xy(1).^2.*xy(2) - xy(2).^3];   
% Jacobian of f
 J = @(xy) [3*xy(1).^2 - 3*xy(2).^2,  -6*xy(1)*xy(2); 
            6*xy(1).*xy(2), 3*xy(1).^2 - 3*xy(2).^2];
        
% Initial guesses
xy = [2; 1];

xy = newton2d(f,J,xy);

x = xy(1); y = xy(2);
disp('parameter values:')

parameter values:

disp(['x = ',num2str(x), '  y = ', num2str(y)])

x = 1  y = -1.0223e-17

disp('function values:' )

function values:

disp(f(xy))

   1.0e-16 *

         0
   -0.3067

function xy = newton2d(f,J,xy)
  
    tol = 10^-6;
    err = 1;
    maxits = 100;
    its = 0;
    while (err>tol) && (its<maxits)
         xyold = xy;
         xy = xyold - J(xyold)\f(xyold);
         err = norm(xy-xyold);
         its = its+1;
    end

end