System of nonlinear equations and Newton's basins of attraction

In [1]:
import numpy as np;
import matplotlib.pyplot as plt;
plt.rcParams.update({"text.usetex":True});
%config InlineBackend.figure_format = "svg"
  • Let us look at following system of non-linear equations
$$x\exp(x+y) - 2 = 0 $$$$xy - 0.1\exp(-y) = 0$$
  • We first solve it using fsolve
In [2]:
import scipy.optimize as sco
In [3]:
def sysfunc(x):
    return [x[0]*np.exp(x[0] + x[1])-2, x[0]*x[1] - 0.1*np.exp(-x[1])];

sol = sco.fsolve(sysfunc, [1, 1]); #provide the system of equations and guess value
print(sol)

[0.80212468 0.11151373]
In [4]:
re1, re2 = sysfunc(sol); 
print(re1, re2)

4.884981308350689e-15 3.788636071533347e-15
  • Let's perform fixed point iterations to solve the same problem
  • Let us re-arrange the first equation as: $$ x = \frac{2}{\exp{(x+y)}}$$ and the second equation as: $$ y = \frac{0.1 \exp{(-y)}}{x}$$
  • We first plot the two functions
    • Since the functions are implicit what we can do is plot contour of function
    • The point where the two contours meet is the root
  • While updating the $x^{(k+1)}$ value we update the value using a weighted sum to dampen out the oscillations
In [5]:
def g(x): # define the system of equations as function
    return np.array([2*np.exp(-(x[0] + x[1])), 0.1/x[0]*np.exp(-x[1])]); #convert list to numpy array

# let's first plot the two implicit functions
x = np.linspace(0, 1);
y = np.linspace(0, 1);
[X, Y] = np.meshgrid(x, y);
F1 = X*np.exp(X+Y) - 2;
F2 = X*Y - 0.1*np.exp(-Y);
plt.contour(X, Y, F1, [0]);# plot the first function
plt.contour(X, Y, F2, [0]);# plot the second function

Niter = 100;
x = np.zeros((Niter, 2));
x[0,:] = np.array([1, 1]); # guess value 1,1
count = 1;

for i in np.arange(1, Niter):
    x[count,:] = 0.6*x[count-1,:]+ 0.4*g(x[count-1,:]);
    count = count + 1;

# plot all the iterated values
plt.plot(x[:,0], x[:,1], '-ok', linewidth=1);

Newton-Raphson method implementation

  • In Newton-Raphson method we have the following formulation: $$ x^{(k+1)} = x^{(k)} - \frac{f^{(k)}}{f`^{(k)}}$$
  • Let's take the case of $$ f(x) = 1-x^2 $$
In [6]:
def f(x): # function
    return 1-x**2;

def dfdx(x): # derivative of function
    return -2*x;

xa = np.linspace(-2, 2);
ya = f(xa);
plt.plot(xa, ya);
plt.plot(xa, 0*xa); # plot the y=0 line
Niter = 10;
x = -0.3; # guess value
for i in np.arange(0, Niter):
    plt.plot([x, x], [0, f(x)], '-k');
    x0 = x;
    x = x - f(x)/dfdx(x);
    plt.plot([x0, x], [f(x0), 0], '-r');

Basins of attraction

  • We shift our focus to roots in the complex plane
  • Let's take $$ f(z) = 1-z^4 $$
  • We define an error threshold as the convergence criteria
    • Error is defined as $$ err(z^{(k+1)}) = \frac{z^{(k+1)} - z^{(k)}}{z^{(k+1)}} $$
  • We also define the maximum number of iterations
  • We define a matrix $\textrm{C}$ that stores the number of iteration required to reach the root for the corresponding guess value in argand plane
  • We then plot the matrix $\textrm{C}$ to observe the regions of convergence
In [7]:
def f(z):
    return 1-z**4;

def dfdz(z):
    return -4*z**3;

def iterf(z):
    Nitermax = 100;
    count = 1;
    err = 1;
    err_threshold = 1e-6;
    while err>err_threshold and count < Nitermax:
        zn = z - f(z)/dfdz(z);
        err = np.abs(1-zn/z);
        z = zn;
        count = count +1;
    return z, count

x = np.linspace(-1.50, 1.5, 200);
y = np.linspace(-1.5, 1.5, 200);
[X, Y] = np.meshgrid(x, y); 
Z = X+1j*Y; # The complex plane

C = np.zeros(np.shape(X)); # stores the number of iteration required to reach the defined convergence criteria of each point in the argand plane
for i in np.arange(0, np.shape(Z)[0]):
    for j in np.arange(0, np.shape(Z)[1]):
        z, C[j,i] = iterf(Z[j,i]);
        #print(Z[j, i], np.round(z, 1), C[j,i] );
        
plt.imshow(C); plt.colorbar();
ax = plt.gca();
ax.set_aspect(1);
  • We observe band of zones where convergence takes more number of iterations than other regions.
  • The dark blue zones are the zones of attraction

Till now we have observed the regions of convergence. Now let's look at which root does the iteration converge to depending upon the initial guess value.

  • For each point in the argand plane we assign it a value corresponding to the root it converges to.
    • We can create a root array and then assign the index of the particular root the point converges to
  • If the root doesn't converge after max number of iterations we assign it a NaNvalue. NaN means not a number.
  • Whenever we encounter a new root, we append it to our root array
In [8]:
def f(z):
    return 1-z**4;

def dfdz(z):
    return -4*z**3;

def iterf(z, Nitermax):
    count = 1;
    err = 1;
    err_threshold = 1e-6;
    while err>err_threshold and count < Nitermax:
        zn = z - f(z)/dfdz(z);
        err = np.abs(1-zn/z);
        z = zn;
        count = count +1;
    return z, count

x = np.linspace(0.25, 0.28, 500);
y = np.linspace(0.25, 0.28, 500);
[X, Y] = np.meshgrid(x, y);
Z = X+1j*Y;

Nitermax = 200;
total_roots = 4;
root_array = np.zeros(total_roots, dtype=complex);
root_count = 0;

C = np.zeros(np.shape(X));
for i in np.arange(0, np.shape(Z)[0]):
    for j in np.arange(0, np.shape(Z)[1]):
        z, count = iterf(Z[j,i], Nitermax);
        z = np.round(z, 1);
        if count < Nitermax: 
            if z not in root_array:
                root_array[root_count] = z;
                root_count = root_count + 1;
        
        if z in root_array:
            C[j,i] = np.where(root_array == z)[0];
        else:
            C[j,i] = np.nan;
In [9]:
plt.imshow(C); plt.colorbar();
In [10]:
print(root_array)

[-1.+0.j -0.-1.j -0.+1.j  1.-0.j]
In [ ]: