Exercise 6.1.14

For a system of equations, e.g. 

$$\begin{array}{llll}\hfill f_1(x,y)& =0& \hfill & \\ \hfill f_2(x,y)& =0& \hfill & \\ \hfill & & \hfill \end{array}$$

If we want to find the roots $(x,y)$ for the two functions, assuming $(x_0,y_0)$ is our initial guess, approximate the equation with Taylor expansion we have

$$\begin{array}{llll}\hfill f_1(x,y)& =f_1(x_0,y_0)+\frac{\partial f_1}{\mathit{\partial x}}{\text{|}}_{(x_0,y_0)}(x-x_0)+\frac{\partial f_1}{\mathit{\partial y}}{\text{|}}_{(x_0,y_0)}(y-y_0)& \hfill & \\ \hfill f_2(x,y)& =f_2(x_0,y_0)+\frac{\partial f_2}{\mathit{\partial x}}{\text{|}}_{(x_0,y_0)}(x-x_0)+\frac{\partial f_2}{\mathit{\partial y}}{\text{|}}_{(x_0,y_0)}(y-y_0)& \hfill & \\ \hfill & & \hfill \end{array}$$

Written in matrix format we have

$$\begin{array}{llll}\hfill \left[\begin{array}{c}\hfill f_1(x,y)\hfill \\ \hfill f_2(x,y)\hfill \end{array}\right]& =\left[\begin{array}{c}\hfill f_1(x_0,y_0)\hfill \\ \hfill f_2(x_0,y_0)\hfill \end{array}\right]+\left[\begin{array}{cc}\hfill \frac{\partial f_1}{\mathit{\partial x}}\hfill & \hfill \frac{\partial f_1}{\mathit{\partial y}}\hfill \\ \hfill \frac{\partial f_2}{\mathit{\partial x}}\hfill & \hfill \frac{\partial f_2}{\mathit{\partial y}}\hfill \end{array}\right]\left(\left[\begin{array}{c}\hfill x\hfill \\ \hfill y\hfill \end{array}\right]-\left[\begin{array}{c}\hfill x_0\hfill \\ \hfill y_0\hfill \end{array}\right]\right)& \hfill & \\ \hfill & & \hfill \end{array}$$

As we are seeking the roots of functions $f_1,f_2$, if $(x,y)$ is the point we are looking for, set $f_1(x,y)=0,f_2(x,y)=0$, we have

$$\begin{array}{llll}\hfill \left[\begin{array}{c}\hfill f_1(x_0,y_0)\hfill \\ \hfill f_2(x_0,y_0)\hfill \end{array}\right]+\left[\begin{array}{cc}\hfill \frac{\partial f_1}{\mathit{\partial x}}\hfill & \hfill \frac{\partial f_1}{\mathit{\partial y}}\hfill \\ \hfill \frac{\partial f_2}{\mathit{\partial x}}\hfill & \hfill \frac{\partial f_2}{\mathit{\partial y}}\hfill \end{array}\right]\left(\left[\begin{array}{c}\hfill x\hfill \\ \hfill y\hfill \end{array}\right]-\left[\begin{array}{c}\hfill x_0\hfill \\ \hfill y_0\hfill \end{array}\right]\right)& =0& \hfill & \\ \hfill \left[\begin{array}{c}\hfill x\hfill \\ \hfill y\hfill \end{array}\right]& =\left[\begin{array}{c}\hfill x_0\hfill \\ \hfill y_0\hfill \end{array}\right]-{\left[\begin{array}{cc}\hfill \frac{\partial f_1}{\mathit{\partial x}}\hfill & \hfill \frac{\partial f_1}{\mathit{\partial y}}\hfill \\ \hfill \frac{\partial f_2}{\mathit{\partial x}}\hfill & \hfill \frac{\partial f_2}{\mathit{\partial y}}\hfill \end{array}\right]}^{-1}\left[\begin{array}{c}\hfill f_1(x_0,y_0)\hfill \\ \hfill f_2(x_0,y_0)\hfill \end{array}\right]& \hfill & \\ \hfill \left[\begin{array}{c}\hfill x\hfill \\ \hfill y\hfill \end{array}\right]& =\left[\begin{array}{c}\hfill x_0\hfill \\ \hfill y_0\hfill \end{array}\right]-J^{-1}\left[\begin{array}{c}\hfill f_1(x_0,y_0)\hfill \\ \hfill f_2(x_0,y_0)\hfill \end{array}\right]& \hfill & \end{array}$$

Apply the formula to the equations in problem 14, we compute the Jacobian matrix first:

$$\begin{array}{llll}\hfill J& =\left[\begin{array}{cc}\hfill \frac{\partial f_1}{\mathit{\partial u}}\hfill & \hfill \frac{\partial f_1}{\mathit{\partial v}}\hfill \\ \hfill \frac{\partial f_2}{\mathit{\partial u}}\hfill & \hfill \frac{\partial f_2}{\mathit{\partial v}}\hfill \end{array}\right]& \hfill & \\ \hfill & =\left[\begin{array}{cc}\hfill 3u^2\hfill & \hfill -1\hfill \\ \hfill -1\hfill & \hfill 3v^2\hfill \end{array}\right]& \hfill & \\ \hfill & & \hfill \end{array}$$

Thus $J^{-1}=\frac{1}{9u^2{v}^2-1}\left[\begin{array}{cc}\hfill 3v^2\hfill & \hfill 1\hfill \\ \hfill 1\hfill & \hfill 3u^2\hfill \end{array}\right]$

import numpy as np 
import math 
 
def newton_m(f, J, x0, maxit=500): 
   x = x0 
   xs, ys = [x0], [f(x0)] 
 
   try: 
       for it in np.arange(maxit): 
           J1 = np.linalg.inv(J(x)) 
           x = x - np.matmul(J1, f(x)) 
           if np.linalg.norm(f(x)) <= 1.0e-12: 
               #print(’Find the root:’, x, f(x), ’ at iteration: ’, it) 
               break 
           xs.append(x) 
           ys.append(f(x)) 
   except Exception as e: 
       return x, np.array([math.inf, math.inf]), None, None, it 
   if it >= maxit: 
       print("Exhausted␣all␣", it, "␣iterations.") 
   return x, f(x), xs, ys, it 
 
def f14(x): 
   u = x[0] 
   v = x[1] 
   res = np.array([u**3 -v, v**3 - u]) 
   return res 
 
def J14(x): 
   u = x[0] 
   v = x[1] 
   res = np.array([[3*u**2, -1], [-1, 3*v**2]]) 
   return res 
 
print("####␣Problem␣VI.1.14") 
solutions = [np.array([0,0]), np.array([1,1]), np.array([-1, -1]), np.array([math.inf, math.inf])] 
data_types = {0:[], 1:[], 2:[], 3:[]} 
xa = np.arange(-1, 1, 0.01) 
ya = np.arange(-1, 1, 0.01) 
#xv, yv = np.meshgrid(xa, ya) 
for i,x in enumerate(xa): 
   for j, y in enumerate(ya): 
       x0 = (x,y) 
       #print("With initial values: ", x0) 
       v, fv, xs, ys, it = newton_m(f14, J14, x0) 
       for idx, sol in enumerate(solutions): 
           if np.isinf(fv[0]) and np.isinf(fv[1]): 
               data_types[3].append(x0) 
           elif np.allclose(v, sol): 
               data_types[idx].append(x0) 
#data_types

\#\#\#\# Problem VI.1.14
    

\ipykernel_launcher.py:26:
RuntimeWarning: overflow encountered in double_scalars
\ipykernel_launcher.py:11:
RuntimeWarning: overflow encountered in matmul
  # This is added back by InteractiveShellApp.init_path()
\ipykernel_launcher.py:11:
RuntimeWarning: invalid value encountered in matmul
  # This is added back by InteractiveShellApp.init_path()
\ipykernel_launcher.py:32:
RuntimeWarning: overflow encountered in double_scalars
\ipykernel_launcher.py:26:
RuntimeWarning: invalid value encountered in double_scalars

import matplotlib.pyplot as plt 
colors = [’y’,’g’,’b’, ’r’] # (0,0), (1,1), (-1,-1), (infty, infty) 
for k, x0s in data_types.items(): 
   #print(x0s) 
   xs = [x for x,y in x0s] 
   ys = [y for x,y in x0s] 
   ys = np.array(ys) 
   plt.scatter(xs, ys, c = colors[k], marker=’.’, s = 0.6)