How do I iterate/increment inorder to find the roots for a system of equations with two variables using the Newton Method

Question

As stated in the title I am interested in finding my roots for a system of two equations and two variables i.e:

$f(x,y)=0$ $g(x,y)=0$

for a the functions:

$f(x,y) =x(1+x^2)-1 $

&

$g(x,y) = y(1+x^2)-2$

How do I iterate properly, since I suspect that my code breaks apart in that regard?

I've confirmed by computing this manually in the command prompt and by the looks of it I get fairly good approximations for the root after 4 iterations. When the tolerance is sufficently close to 0.

Other than that I've just changed vector lengths and various values. Nothing worth mentioning.

function [x,y] = NewtonMeth(f, g, a, b)

n=100; %arbitrary amount of iterations
x=zeros(1,n+1); %initializing my vector for x values 
y=zeros(1,n+1);  %initializing my vector for y values
y(1)=b; %setting my first element as the intial guess
x(1)=a;


tol = 0.00001; %tolerance, denoting when to stop my for loop.

for i=1:n
%utilizing a previous function for partial derivative:

 f1=PD(f, x(i), y(i), 1); %partial derivatives for my two functions f & g 
 f2=PD(f, x(i), y(i), 2);
 g1=PD(g, x(i), y(i), 1);
 g2=PD(g, x(i), y(i), 2);

 j1= [f(x(i),y(i)) f2; g(x(i),y(i)) g2]; %separating my jacobian into amatrix containing the relevant elements
 j2 = [f1 f2; g1 g2];%second jacobian matrix
 j = det(j1)./det(j2); %here I evaluate the jacobian, with accordance to the formula for the newton method for system of equations 

   %with two variables. 




   if abs(f(x(i),y(i)))>=tol %Said tolerance which should limit the amount of iterations


      x(i+1)= x(i)-j; %incrementing the x values inorder to approximate the root 

      y(i+1)=y(i)-j; %incrementing the y values inorder to approximate the root



     end
   end
end

Here is my function for taking the partial derivative:

function derivative = PD(f, a, b, i)

h = 0.000001;
fn=zeros(1,2);

if i == 1
fn(i) = (f(a+h,b)-f(a,b))/h;


 elseif i==2

    fn(i) = (f(a,b+h)-f(a,b))/h;

  end

 derivative = fn(i);

end

Choosing the intial guess for a=0.2 and b= 1.8 (as it's done in my calc 3 book)

I should get:

|x(1) = 0.200000|y(1)=1.800000|f(x(1),y(1))=-0.152000|g(x(1),y(1))=-0.128000|
|x(2) = 0.216941|y(2)=1.911349|f(x(2),y(2))= 0.009481|g(x(2),y(2))= 0.001303|
|x(3) = 0.214827|y(3)=1.911779|f(x(3),y(3))=-0.000003|g(x(3),y(3))= 0.000008|
|x(4) = 0.214829|y(4)=1.911769|f(x(4),y(4))= 0.000000|g(x(4),y(4))= 0.000000|

Answer 1

Without considering the coding aspects, solving for $2$ variables is basically the same as with one variable except that the expansions are to be done for $x$ and $y$.

So, locally, we have $$f(x,y)=f(a,b)+f_x(a,b)(x-a)+f_y(a,b)(y-b)=0$$ $$g(x,y)=g(a,b)+g_x(a,b)(x-a)+g_y(a,b)(y-b)=0$$ and you solve two linear equations in $(x-a)$ and $(y-b)$. When solved, replace $a$ and $b$ by the new values of $x$ and $y$.

You compute the derivatives according to $$f_x(a,b)=\frac{f(a+\Delta a,b)-f(a,b)}{\Delta a}\qquad f_y(a,b)=\frac{f(a,b+\Delta b)-f(a,b)}{\Delta b}$$ $$g_x(a,b)=\frac{g(a+\Delta a,b)-g(a,b)}{\Delta a}\qquad g_y(a,b)=\frac{g(a,b+\Delta b)-g(a,b)}{\Delta b}$$

If you want, for sure, you can make $\Delta a=\Delta b=h$.

In any manner, I think that yo have typo's somewhere in $f(x,y)$ since, as written, $f(0.2,1.8)=-0.792$ and not $-0.152$.