import numpy as np;
import matplotlib.pyplot as plt;
plt.rcParams.update({"text.usetex":True}); # change inline style to latex
%config InlineBackend.figure_format = "svg" # change figure type to svg
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.
x = np.linspace(-2,5,200); # create the points between which the root might exist
y = np.exp(x)-3*x**2;
plt.plot(x, y, label="$f(x) = \exp(x)-3x^2$"); plt.plot(x, 0*x);
plt.ylim([-2, 2]);
plt.legend(fontsize=16);
plt.xlabel("$x$", fontsize=16);
plt.ylabel("$y = f(x)$", fontsize=16)
plt.title("Intersection points of $f(x)$ with X-axis are the roots");
x0 = 2.0;
Niter = 20;
x = np.zeros(Niter); # Initialization. Stores the value of ith iteration
x[0] = x0;
print("Iteration: %d \t x: %f"%(0,x[0])); # This is for printing the 0th iterate - the intial value
# This function returns the RHS of the equation
def fx(x):
return np.log(3) + 2*np.log(x);
# This is the loop where the values are updated
# We will print out the iteration number and updated value
for i in np.arange(1,Niter):
x[i] = fx(x[i-1]);
print("Iteration: %d \t x: %3.16f"%(i,x[i]));
# Print the array which contains all the elements
print(x);
plt.plot(np.arange(0,Niter), x, '-ok');
Let us convert that into a function to plot all the family of curves.
fx(x) takes as an input the value $x$ and returns $\ln 3 + 2\ln x$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.x_0 to the function and plot the different arrays x. def fx(x):
return np.log(3) + 2*np.log(x);
def perform_iterations(x0, Niter):
x = np.zeros(Niter); # Initialization
x[0] = x0;
for i in np.arange(1,Niter):
x[i] = fx(x[i-1]);
return x
Niter = 20;
x0_a = np.linspace(1,5, 20);
for x0 in x0_a:
x = perform_iterations(x0, Niter);
plt.plot(np.arange(0, Niter), x);
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
x0 = 2.0;
Niter = 20;
x = np.zeros(Niter); # Initialization
x[0] = x0;
def fx(x):
return 1/3**(0.5)*np.exp(x/2);
for i in np.arange(1,Niter):
x[i] = fx(x[i-1]);
print(x);
We find that this representation, too, converges to a root. In the same way, the third root can also be found by starting with an appropriate initial guess. Try it out yourself!
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$$
x = np.linspace(-2, 4);
y = np.exp(-x) - 3*x;
plt.plot(x,y); plt.plot(x, 0*x);
plt.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 .
x0 = 2.0; # Guess value
Niter = 20; # Number of iterations
x = np.zeros(Niter); # Initialization
x[0] = x0; # Setting the initial element of the array x
def fx(x):
return 1/3*np.exp(-x);
for i in np.arange(1,Niter):
x[i] = fx(x[i-1]);
print(x, fx(x[Niter-1]));
print(np.isclose(x, fx(x[Niter-1])));
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.
x0 = 2.0;
Niter = 20;
x = np.zeros(Niter); # Initialization
xx = np.zeros(Niter);
x[0] = x0;
xx[0] = x0;
def fx(x): # THis function diverges
return np.exp(-x)-2*x;
def gx(x):
return 1/3*np.exp(-x); # THis function converges
for i in np.arange(1,Niter):
x[i] = fx(x[i-1]);
xx[i] = gx(xx[i-1]);
f(x) divergesplt.plot(np.arange(0, Niter), x, np.arange(0, Niter), xx)
plt.ylim([-2, 2])
As seen from the curve of x vs Niter, the orange 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:
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.
Niter = 10;
x = np.zeros(Niter);
xg = 1.1;
def fx(x):
return np.log(3) + 2*np.log(x);
xp = np.linspace(0.2, 5, 100);
yp = fx(xp);
plt.plot(xp, yp, label='$\\ln 3 + 2\\ln x$');
plt.plot(xp, xp, label='$y = x$');
plt.legend(fontsize=16);
plt.xlabel("$x$", fontsize=16);
plt.ylabel("$x, f(x)$", fontsize=16);
plt.xticks(fontsize=16); plt.yticks(fontsize=16);
ax = plt.gca(); ax.set_aspect(0.6);
for i in np.arange(1, Niter):
xn = fx(xg);
plt.plot(xg, fx(xg), 'ok', markerfacecolor='black', markersize=2);
x2 = xn; y2 = x2;
plt.plot(x2, y2, 'ok', markerfacecolor='black', markersize=2 );
plt.arrow(xg, fx(xg), x2-xg, y2-fx(xg), length_includes_head=True, linewidth=0.7, color='b');
plt.arrow(x2, y2, x2-x2, fx(x2)-y2, length_includes_head=True, linewidth=0.7, color='r');
xg = xn;
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}{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.
Niter = 10;
x = np.zeros(Niter);
xg =3.5;
def fx(x):
return 1/3**(0.5)*np.exp(x/2);
xp = np.linspace(0.2, 5, 100);
yp = fx(xp);
plt.plot(xp, yp, label='$\\frac{1}{3}\\exp\left(\\frac{x}{2}\\right)$');
plt.plot(xp, xp, label='$y = x$');
plt.legend(fontsize=16);
plt.xlabel("$x$", fontsize=16);
plt.ylabel("$x, f(x)$", fontsize=16);
plt.xticks(fontsize=16); plt.yticks(fontsize=16);
ax = plt.gca(); ax.set_aspect(0.6);
for i in np.arange(1, Niter):
xn = fx(xg);
plt.plot(xg, fx(xg), 'ok', markerfacecolor='black', markersize=2);
x2 = xn; y2 = x2;
plt.plot(x2, y2, 'ok', markerfacecolor='black', markersize=2 );
plt.arrow(xg, fx(xg), x2-xg, y2-fx(xg), length_includes_head=True, linewidth=0.7, color='b');
plt.arrow(x2, y2, x2-x2, fx(x2)-y2, length_includes_head=True, linewidth=0.7, color='r');
xg = xn;
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$$.
Niter = 10;
x = np.zeros(Niter);
xg =0.25;
def fx(x):
return np.exp(-x)-2*x;
xp = np.linspace(-5, 5, 100);
yp = fx(xp);
plt.plot(xp, yp, label='RHS function');
plt.plot(xp, xp, label='$y = x$');
plt.xlim(0, 0.5)
plt.ylim(-0.5,2)
for i in np.arange(1, Niter):
xn = fx(xg);
plt.plot(xg, fx(xg), 'ok', markerfacecolor='black', markersize=2);
x2 = xn; y2 = x2;
plt.plot(x2, y2, 'ok', markerfacecolor='black', markersize=2 );
plt.arrow(xg, fx(xg), x2-xg, y2-fx(xg), length_includes_head=True, linewidth=0.7, color='b');
plt.arrow(x2, y2, x2-x2, fx(x2)-y2, length_includes_head=True, linewidth=0.7, color='r');
xg = xn;
There are a few observables from the plot:
We can use one of several functions available in the submodule scipy.optimize for solving nonlinear equations. Let us demonstrate this below:
import scipy.optimize as sco;
Let us solve the nonlinear equation $\exp(x) - 3x^2 =0$ using the
fsolve function of pythonbrentq function of python; this is a bracketing methodnewton 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 methodIn all the cases the functions are passed to the solver by means of a function handle. This is illustrated below:
def fx(x):
return np.exp(x)-3*x**2;
def fpx(x):
return np.exp(x)-6*x;
sol = sco.fsolve(fx, -2);
print(sol)
sol2 = sco.brentq(fx, -2, 0);
print(sol2)
sol3 = sco.newton(fx, -2, fpx); # Jacobian
print(sol3)