How to use this Matlab code to write another, about Newton's method of system of equations [closed]

Question

Here is a code for applying Newton's for finding roots for $x^2-1$. However, how to use this code or just make some changes and define something more,(like define the column vector of $(x_1,x_2)$) on the code below to find roots for this system of equations: $f(x_1,x_2)=x_1^2+x_2^2-1=0$, $f(x_1,x_2)=5x_1^2-x_2^2-2=0$

also using the newton's method. The formula will be: $x_{k+1}=x_k-(▽f(x_k))^{(-T)}f(x_k)$, where $(▽f(x_k))^{T}$ is the Jacobian matrix. What the code will be?

Here is the code for $x^2-1$ with first trial $x(1)=10$:

$function[]=newton_eg1 %% demonstration of Newton's method %% the derivative of f is nonzero at the solution. x(1)=10; epsilon=1.e-10; MaxIter=100; k=1; %% f=x^2-1 while abs(x(k)^2-1)>epsilon & k <MaxIter x(k+1)=x(k) - (x(k)^2-1)/(2*x(k)); k=k+1; end %%%%%%output solution and error: fprintf('step solution error \n'); for i=1:k err(i)=abs(x(i)-x(k)); fprintf('%i: %e %e\n', i, x(i), err(i)); end fprintf('\n convergence rate: \n'); for i=1:k-3 order(i)=log(err(i+2)/err(i+1)) / log( err(i+1)/err(i) ); fprintf('%i: %e \n', i, order(i)); end $

Answer 1

I wrote an example code for the system you gave. I tested this with Octave, but I don't think Matlab should complain. I hope this helps.

function newton_system()
        eps = 1e-10; itmax = 100;
        X = zeros(2,100); % Store all iteration
        X(:,1) = rand(2,1); % Initialize first step 
        for i=2:itmax
                %x_next = x_old - f'^-1 *f
                %      <=>
                % x_next-x_old =- f'^-1 *f
                %      <=>
                % f'(x_next-x_old) = -f
                %      <=>
                % Ax = b
                %      <=>
                %  x = A\b
                %      <=>
                %   x_next = f'\f + x_old
                X(:,i) = df(X(:,i-1))\-f(X(:,i-1)) + X(:,i-1);
                if abs(f(X(:,i)))<eps
                        X = X(:,1:i); % cut remaining zeros
                        disp("The solution is:")
                        disp(X(:,i));
                        break
                end

        end
        F = f(X)
        ERR2NORM = sqrt((F(1,:))-(F(1,end)).^2+(F(2,:)-F(2,end)).^2);
        semilogy(ERR2NORM)

end

function Y=f(X)
        % note that the (1,:) or (2,:) could all be replaced 
        % by (1) and (2) if it wasn't for the fact that I want to use
        % f(X) in line 26 for all elements stored in the vector
        Y = zeros(size(X));
        Y(1,:) = X(1,:).*X(1,:)+X(2,:).*X(2,:)-1;
        Y(2,:) =5*X(1,:).*X(1,:)-X(2,:).*X(2,:)-2;
end

function dY=df(X)
        dY = zeros(2,2);
        dY(1,1) = 2*X(1); dY(1,2) = 2*X(2);
        dY(2,1) = 10*X(1);dY(2,2) = -2*X(2);

end