Fixed point iterations and cobwebs

In this lecture, we will be looking at solving nonlinear equations. In particular, we will be looking at a very simple idea - fixed point iterations. While fixed point iterations have far reaching consequences for understanding various aspects of nonlinear dynamics as well, it lends an insight into the nature of convergence of roots and the dependency of the convergence behaviour depending on the functional choice.

In [1]:
set(0, "defaultlinelinewidth", 5);
set (0, "defaulttextfontname", "TimesNewRoman")
set (0, "defaulttextfontsize", 20)
set (0, "DefaultAxesFontName", "TimesNewRoman")
set(0, 'DefaultAxesFontSize', 20)

Fixed point iterations in 1D

Let us consider the nonlinear equation $$\exp(x) - 3x^2 = 0, $$ whose roots we are interested it. In order to apply fixed point iterations, we will first cast it in the form of $x = f(x)$. Towards this, we can have multiple possibilities to do so. Let's consider the first possibility as shown below.

We rearrange the equation to write $$x = \ln 3 + 2 \ln x.$$ In case we had the root, the left hand side and right hand side of the equations would match exactly. However, if we guess a value different from the root, then the left hand side and right hand side would not match. In that case we can substitute the guess value in the right hand side function; in this case $\ln 3 + 2\ln x$. This value can then be assigned to $x$. This becomes the updated value of $x$.

We can program the aforementioned discussion in the following form: $$x^{(k+1)} = f(x^{(k)}) $$ where $k$ starts from 0 to $n$, where $n$ is the specified number of iterations. When $k =0$, we are at the first guess point, $x^{(0)}$. There is no universal form for the choice of neither the initial guess nor the number of iterations.

Let's first assess how the plot looks like and where to expect the roots to exist.

In [2]:
x = linspace(-2,5,200);
y = exp(x)-3*x.^2;
#plot(x, y,"linewidth",5, x,0.*x,"linewidth",5);
plot(x, y, x,0.*x);
ylim([-2, 2]); 
h = legend ("exp(x) - 3x^2");
legend boxoff
legend (h, "location", "northwest");
set (h, "fontsize", 20,"fontname", "TimesNewRoman");
set(gca,'FontSize',20,"fontname", "TimesNewRoman");
xlabel("x", "fontsize",20,"fontname", "TimesNewRoman");
ylabel("y = f(x)","fontsize",20,"fontname", "TimesNewRoman"); 
title("Intersection points of f(x) with X-axis are the roots", "fontsize",20,"fontname", "TimesNewRoman");
In [3]:
x0 = 2.0;
Niter = 20;
x = zeros(1,Niter); # Initialization
x(1) = x0;

printf("Iteration: %d \t x: %f\n",0,x(1)); # This is for printing the 0th iterate - the intial value

# This function returns the RHS of the equation
function y = fx(x)
    y = log(3) + 2*log(x);
end
# This is the loop where the values are updated
# We will print out the iteration number and updated value
for i=2:Niter+1
    x(i) = fx(x(i-1));
    printf("Iteration: %d \t x: %3.16f\n",i-1,x(i));
end  
    
# Print the array which contains all the elements
disp(x);

Iteration: 0 	 x: 2.000000
Iteration: 1 	 x: 2.4849066497880004
Iteration: 2 	 x: 2.9190824753987616
Iteration: 3 	 x: 3.2411509809654229
Iteration: 4 	 x: 3.4504693041876902
Iteration: 5 	 x: 3.5756327926507705
Iteration: 6 	 x: 3.6468966195075572
Iteration: 7 	 x: 3.6863654180751197
Iteration: 8 	 x: 3.7078942708438589
Iteration: 9 	 x: 3.7195405556699663
Iteration: 10 	 x: 3.7258125969024354
Iteration: 11 	 x: 3.7291822383889315
Iteration: 12 	 x: 3.7309902298311357
Iteration: 13 	 x: 3.7319596398469521
Iteration: 14 	 x: 3.7324792252860597
Iteration: 15 	 x: 3.7327576577037318
Iteration: 16 	 x: 3.7329068464955490
Iteration: 17 	 x: 3.7329867797911067
Iteration: 18 	 x: 3.7330296056332966
Iteration: 19 	 x: 3.7330525500470131
Iteration: 20 	 x: 3.7330648426595179
 Columns 1 through 8:

   2.0000   2.4849   2.9191   3.2412   3.4505   3.5756   3.6469   3.6864

 Columns 9 through 16:

   3.7079   3.7195   3.7258   3.7292   3.7310   3.7320   3.7325   3.7328

 Columns 17 through 21:

   3.7329   3.7330   3.7330   3.7331   3.7331
In [4]:
plot(0:Niter, x, '-ok');

It can be clearly seen from the plot above that we start from an initial value of 2 at the zeroth iteration. The values are successively updated and eventually after 15 iterations we converge to the root of the equation. The iterated values of $x$ are also shown above.

Wrapping the above code to see the influence of the guess value

Let us convert that into a function to plot all the family of curves.

  • The function fx(x) takes as an input the value $x$ and returns $\ln 3 + 2\ln x$
  • The function perform_iterations takes as an input the initial guess x_0 and the number of iterations Niter. Based on these values the function will perform the iterations and return the array x which contains all the iterated values.
  • In the end we will pass different values of x_0 to the function and plot the different arrays x.
In [5]:
function y = fx(x)
    y = log(3) + 2*log(x);
end

function x = perform_iterations(x0, Niter)
    x = zeros(1,Niter); # Initialization
    x(1) = x0;
    for i=2:Niter+1
        x(i) = fx(x(i-1));
    end
end

Niter = 20;
x0_a = linspace(1,5, 20);

for i=1:length(x0_a)
    x = perform_iterations(x0_a(i), Niter);
    plot(0:Niter, x, "linewidth",5);
    hold on;
    drawnow();
end
set(gca,'FontSize',20,"fontname", "TimesNewRoman");
xlabel("x", "fontsize",20,"fontname", "TimesNewRoman");
ylabel("y = f(x)","fontsize",20,"fontname", "TimesNewRoman");

Another rearrangement of the same nonlinear equation

Let us consider the same nonlinear equation $$\exp(x) - 3x^2 = 0$$ but now we will rearrange as $$x = \frac{1}{\sqrt{3}}\exp\frac{x}{2}$$. Let us reuse the same code above

In [6]:
x0 = 2.0;
Niter = 20;
x = zeros(1,Niter); # Initialization
x(1) = x0;

function y = fx(x)
    y = 1/3^(0.5)*exp(x/2);
end

for i=2:Niter+1
    x(i) = fx(x(i-1));
end
disp(x);

 Columns 1 through 8:

   2.00000   1.56940   1.26541   1.08697   0.99420   0.94913   0.92798   0.91822

 Columns 9 through 16:

   0.91375   0.91171   0.91078   0.91036   0.91017   0.91008   0.91004   0.91002

 Columns 17 through 21:

   0.91001   0.91001   0.91001   0.91001   0.91001

We find that this representation, too, convergence to the root. In the same way, the third root can also be found by starting with an appropriate initial guess. Try it out yourself!

A case where convergence may not occur.

While the previous nonlinear equation and its rearrangement to a fixed point iteration showed promise, it may not always be the case. Consider the nonlinear equation $$\exp(-x) - 3x = 0$$. Before attempting to find out the roots of this equation, let us first look at the curve represented by the equation itself.

In [7]:
x = linspace(-2, 4);
y = exp(-x) - 3*x;
plot(x,y, x, 0.*x);
ylim([-2, 2]);

It is clear from the curve above that there does exist a root somewhere between 0 and 1. Somewhere close to 0.25. So our initial guess for the fixed point iterations can be 0.25. Great!

Let us set up the equation in the form $x = f(x)$. We can rearrange $$\exp(-x) - 3x = 0\qquad \textrm{as}\qquad x = \frac{1}{3}\exp(-x)$$. Let's see what that gives us ...

In [8]:
x0 = 2.0;  # Guess value
Niter = 20; # Number of iterations
x = zeros(1,Niter); # Initialization
x(1) = x0; # Setting the initial element of the array x

function y = fx(x)
    y = 1/3*exp(-x);
end

for i=1:Niter
    x(i+1) = fx(x(i));
end

printf("%g %g",x,fx(x(Niter)));
lia  = ismembertol(x, 0.*x + fx(x(Niter)));

2 0.04511180.31863 0.2423810.261586 0.256610.25789 0.257560.257645 0.2576230.257629 0.2576270.257628 0.2576280.257628 0.2576280.257628 0.2576280.257628 0.2576280.257628 0.257628'perl' is not recognized as an internal or external command,
operable program or batch file.
error: 'ismembertol' undefined near line 1 column 8

Well, the solution does converge to the value of 0.2576... with the choice of casting the iterations. The speed with which convergence is acheived can be probed with the help of the function np.isclose(x, fx(x[last])) which tells us which elements of x are close to the fx(x[last]). After the first 10 odd terms, we see that we do not have much change in the iterated value of x.

Is convergence always guaranteed?

The answer is no. This will be seen from the next case study. In this particular case, we will recast the iterations in the following form: $$\exp(-x) - 3x = 0\qquad \textrm{as}\qquad x = \exp(-x) - 2x$$. Let's investigate both the representations together:

In [9]:
x0 = 2.0;
Niter = 20;
x = zeros(1,Niter); # Initialization
xx = zeros(1,Niter);
x(1) = x0;
xx(1) = x0;

function y= fx(x) # THis function diverges
    y = exp(-x)-2*x;
end

function y= gx(x)
    y =  1/3*exp(-x); # THis function converges
end
for i =1:Niter
    x(i+1) = fx(x(i));
    xx(i+1) = gx(xx(i));
end
In [10]:
plot(0:Niter, x, '-b',0:Niter, xx,'-r')
ylim([-2, 2])

As seen from the curve of x vs Niter, the red curve iterates while the blue curve does not converge at all. It should be clear that the choice of the representation leads to the following two observations:

  • The root to which the iterations converge depend on the form of the RHS cast.
  • The iterations may or may not coverge depending on the form of the RHS cast.

Let us investiage this phenomenon further by explicitly plotting how the iterations proceed to actually modify the guess value to the next value. For this, let us consider first the case that showed convergence - $$x = \ln 3 + 2\ln x$$.

We will plot the $f(x)$, which is the RHS of the above iteration scheme. Apart from this we will also plot $x$. The iterations can be logically split into the two steps: $$y_1 = \ln 3 + 2\ln x$$ followed by $$x_1 = y_1$$ hence the two curves give us the way the plots converge.

In [11]:
Niter = 10;
x = zeros(1,Niter);
xg = 1.1;

function y= fx(x)
    y= log(3) + 2*log(x);
end

xp = linspace(0.2, 5, 100);
yp = fx(xp);
plot(xp, yp, "-;\ln 3 + 2\ln x;");     #plot(x, y, sprintf("-;J_{%d};",i), "linewidth",5);
hold on;
plot(xp, xp, "-;y = x;");
h = legend();
legend (h, "location", "southwest");
set (h, "fontsize", 16,"fontname", "TimesNewRoman");
legend boxoff
xlabel("x");
ylabel("x, f(x)");
#ax = gca(); set_aspect(0.6);

for i=1:Niter
    xn = fx(xg);
    plot(xg, fx(xg), 'ok', markersize=2);
    x2 = xn; y2 = x2;
    plot(x2, y2, 'ok', markersize=2 );
    plot([xg, x2], [fx(xg), y2],'Color', 'blue' )
    plot([x2, x2], [y2, fx(x2)], 'Color','red' )
    xg = xn;
    drawnow();
end

In the above curve, we have started with an initial guess near 1. This eventually iterates towards the root in the vicinity of 3.7. The intersection of the two curves are the roots since at those points of intersection $y = x$ and $y = f(x)$ $\implies$ $x = f(x)$.

Interestingly, despite starting the guess near the root near 1, we have coverged to the other root. We did, however, see that the second representation of the nonlinear equation, i.e. $$ \frac{1}{\sqrt{3}}\exp\left(\frac{x}{2}\right) $$ does give us the other root. Let's investigate the sequence of convergence for this particular case. It is easy for us to simply reuse the above code block and change the function to the desired one.

In [12]:
Niter = 10;
x = zeros(1,Niter);
xg = 3.5;

function y= fx(x)
    y= 1/3^(0.5)*exp(x/2);
end

xp = linspace(0.2, 5, 100);
yp = fx(xp);
plot(xp, yp, "-;\ln 3 + 2\ln x;");     #plot(x, y, sprintf("-;J_{%d};",i), "linewidth",5);
hold on;
plot(xp, xp, "-;y = x;");
h = legend();
legend (h, "location", "northwest");
set (h, "fontsize", 16,"fontname", "TimesNewRoman");
legend boxoff
xlabel("x");
ylabel("x, f(x)");
#ax = gca(); set_aspect(0.6);

for i=1:Niter
    xn = fx(xg);
    plot(xg, fx(xg), 'ok', markersize=2);
    x2 = xn; y2 = x2;
    plot(x2, y2, 'ok', markersize=2 );
    plot([xg, x2], [fx(xg), y2],'Color', 'blue' )
    plot([x2, x2], [y2, fx(x2)], 'Color','red' )
    xg = xn;
    drawnow();
end

This particular sequence does depict convergence towards the other root. This happens despite choosing an initial guess quite close to the other root. If we compare the two curves, it is clear that in one case the curvature is concave down while in the other it is concave up. Regardless of that, the important thing is the slope of the curve at the root. In the slope of $f(x)$ at the root $X$ is such that $|df/dx|_X|<1$, then the iterations converge towards that root. This explains why deciding the shape of $f(x)$ has such a critical influence on the direction of convergence.

With this in mind, let's investigate the case where the iterations were seen to diverge despite choosing an initial guess quite close to the root. To recall, the iterations were setup in the following manner: $$x = \exp(-x)-2x$$.

In [13]:
Niter = 10;
x = zeros(1,Niter);
xg =0.25;

function y= fx(x)
    y = exp(-x)-2*x;
end

xp = linspace(-5, 5, 100);
yp = fx(xp);
plot(xp, yp, "-;RHS function;");
hold on;
plot(xp, xp, "-; y = x;");
xlim([0, 0.5]);
ylim([-0.5,2]);

for i=1:Niter
    xn = fx(xg);
    plot(xg, fx(xg), 'ok', markersize=2);
    x2 = xn; y2 = x2;
    plot(x2, y2, 'ok', markersize=2 );
    plot([xg, x2], [fx(xg), y2],'Color', 'blue' )
    plot([x2, x2], [y2, fx(x2)], 'Color','red' )
    xg = xn;
    drawnow();
end

There are a few observables from the plot:

  • We have initialized really close to the root, despite that the root seems to fly away.
  • The negative slope of the curve leads to a spiral nature of iterations.
  • The magnitude of the slope at the root is larger than 1 at the root (verify this yourself). This causes the instability.
  • The figure looks like a spider cobweb!

Using Python inbuilt libraries

We can use one of several functions available in the submodule scipy.optimize for solving nonlinear equations. Let us demonstrate this below:

Let us solve the nonlinear equation $\exp(x) - 3x^2 =0$ using the

  • fsolve function of python
  • brentq function of python; this is a bracketing method
  • newton function of python; this is the classic Newton-Raphson method. We must supply a Jacobian for this. This is essentially the derivative (analytical) of the function. In case we do not specify the derivative, the secant method is employed in place of the Newton-Raphson method

In all the cases the functions are passed to the solver by means of a function handle. This is illustrated below:

In [14]:
function y= fx(x)
    y= exp(x)-3*x.^2; 
end

function y= fpx(x)
    y= exp(x)-6*x; 
end

[sol, res] = fsolve (@fx, -2);
disp(sol)

sol2 = fzero(@fx, -2);
disp(sol2)

-0.45896
-0.45896