System of nonlinear equations and Newton's basins of attraction

In [1]:
# Set default parameters for the plots
set(0, "defaultlinelinewidth", 5);
set (0, "defaulttextfontname", "TimesNewRoman")
set (0, "defaulttextfontsize", 20)
set (0, "DefaultAxesFontName", "TimesNewRoman")
set(0, 'DefaultAxesFontSize', 20)
  • 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]:
function y = sysfunc(x)
    y(1) = x(1)*exp(x(1) + x(2))-2;
    y(2) = x(1)*x(2) - 0.1*exp(-x(2));
end
[sol, res] = fsolve (@sysfunc, [1; 1]);
disp(sol)
disp(res)

   0.80212
   0.11151
  -0.000000000015377  -0.000000000080634
In [3]:
res = sysfunc(sol); 
printf("%g \t %g",res(1), res(2))

-1.53766e-11 	 -8.06345e-11
  • 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 [4]:
function y= g(x)
    y = [2*exp(-(x(1) + x(2))), 0.1/x(1)*exp(-x(2))];
end

x = linspace(0, 1);
y = linspace(0, 1);
[X, Y] = meshgrid(x, y);
F1 = X.*exp(X+Y) - 2;
F2 = X.*Y - 0.1.*exp(-Y);
contour(X, Y, F1, [0, 0]);
hold on;
contour(X, Y, F2, [0, 0]);

Niter = 100;
x = zeros(Niter,2);
x(1,:) = [1, 1];

for i=2:Niter
    x(i,:) = 0.6*x(i-1,:)+ 0.4*g(x(i-1,:));
end
plot(x(:,1), x(:,2), '-ok', 'markerfacecolor', 'black');
hold off;

disp(x(end,end))

 0.11151

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 [5]:
function y= f(x)
    y = 1-x.^2;
end

function dy= dfdx(x)
    dy= -2*x;
end

xa = linspace(-2, 2);
ya = f(xa);
plot(xa, ya);
hold on;
plot(xa, 0.*xa);
Niter = 10;
x = -0.3;
for i=0:Niter
    plot([x, x], [0, f(x)], '-k');
    x0 = x;
    x = x - f(x)/dfdx(x);
    plot([x0, x], [f(x0), 0], '-r');
end

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 [6]:
function y= f(z)
    y =  1-z.^4;
end

function dz = dfdz(z)
    dz =  -4*z.^3;
end

function [z,count] =  iterf(z1)
    Nitermax = 100;
    count = 1;
    err = 1;
    err_threshold = 1e-6;
    while (err>err_threshold) && (count < Nitermax)
        zn = z1 - f(z1)/dfdz(z1);
        err = abs(1-zn/z1);
        z1 = zn;
        count = count +1;
    end
    z = z1;
end

x = linspace(-1.5, 1.5, 200);
y = linspace(-1.5, 1.5, 200);
[X, Y] = meshgrid(x, y);
Z = X+1j*Y;

C = zeros(size(X));
size(Z)(1)

ans =  200
In [7]:
for i=1:size(Z)(1)
    for j=1:size(Z)(2)
        [z, C(j,i)] = iterf(Z(j,i));
    end
end
In [8]:
image(C); colorbar;
  • 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 [9]:
function y= f(z)
    y =  1-z.^4;
end

function dz = dfdz(z)
    dz =  -4*z.^3;
end

function [z,count] =  iterf(z1, Nitermax)
    count = 1;
    err = 1;
    err_threshold = 1e-6;
    while (err>err_threshold) && (count < Nitermax)
        zn = z1 - f(z1)/dfdz(z1);
        err = abs(1-zn/z1);
        z1 = zn;
        count = count +1;
    end
    z = z1;
end


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

Nitermax = 200;
total_roots = 4;
root_array = zeros(total_roots,1);
root_count = 1;

C = zeros(size(X));

for i=1:size(Z)(1)
    for j=1:size(Z)(2)
        [z, count] = iterf(Z(j,i), Nitermax);
        z = round(z);
        if count < Nitermax 
            if ~ismember(1, (z==root_array))
                root_array(root_count) = z;
                root_count = root_count + 1;
            end
        end
        if ismember(1, (z== root_array))
            C(j,i) = find(1== (z==root_array));
        else
            C(j,i) = NaN;
        end
    end
end
In [12]:
pcolor(flipud(C)); colorbar(); shading interp;
  • We make use of flipud command because in image the top left corner is (0,0)
In [11]:
root_array

root_array =

   0 + 1i
   0 - 1i
  -1 + 0i
   1 + 0i
In [ ]: