In [1]:
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
%matplotlib inline
from ipywidgets import interact
from ipywidgets import widgets
sym.init_printing()
from scipy import optimize
Hello again! In this document we're going to learn how to find a 1D equation's solution using numerical methods. First, let's start with the definition of a root:
Definition: The function $f(x)$ has a root in $x = r$ if $f(r) = 0$.
An example: Let's say we want to solve the equation $x + \log(x) = 3$. We can rearrange the equation: $x + \log(x) - 3 = 0$. That way, to find its solution we can find the root of $f(x) = x + \log(x) - 3$. Now let's study some numerical methods to solve these kinds of problems.
Defining a function $f(x)$
In [2]:
f = lambda x: x+np.log(x)-3
Finding $r$ using sympy
In [3]:
y = sym.Symbol('y')
fsym = lambda y: y+sym.log(y)-3
r_all=sym.solve(sym.Eq(fsym(y), 0), y)
r=r_all[0].evalf()
print(r)
print(r_all)
In [4]:
def find_root_manually(r=2.0):
x = np.linspace(1,3,1000)
plt.figure(figsize=(8,8))
plt.plot(x,f(x),'b-')
plt.grid()
plt.ylabel('$f(x)$',fontsize=16)
plt.xlabel('$x$',fontsize=16)
plt.title('What is r such that $f(r)='+str(f(r))+'$? $r='+str(r)+'$',fontsize=16)
plt.plot(r,f(r),'k.',markersize=20)
plt.show()
interact(find_root_manually,r=(1e-5,3,1e-3))
Out[4]:
The bisection method finds the root of a function $f$, where $f$ is a continuous function. If we want to know if this has a root, we have to check if there is an interval $[a,b]$ for which $f(a)\cdot f(b) < 0$. When these 2 conditions are satisfied, it means that there is a value $r$, between $a$ and $b$, for which $f(r) = 0$. To summarize how this method works, start with the aforementioned interval (checking that there's a root in it), and split it into two smaller intervals $[a,c]$ and $[c,b]$. Then, check which of the two intervals contains a root. Keep splitting each "eligible" interval until the algorithm converges or the tolerance is surpassed.
In [5]:
def bisect(f, a, b, tol=1e-8):
fa = f(a)
fb = f(b)
i = 0
# Just checking if the sign is not negative => not root necessarily
if np.sign(f(a)*f(b)) >= 0:
print('f(a)f(b)<0 not satisfied!')
return None
#Printing the evolution of the computation of the root
print(' i | a | c | b | fa | fc | fb | b-a')
print('----------------------------------------------------------------------------------------')
while(b-a)/2 > tol:
c = (a+b)/2.
fc = f(c)
print('%2d | %.7f | %.7f | %.7f | %.7f | %.7f | %.7f | %.7f' %
(i+1, a, c, b, fa, fc, fb, b-a))
# Did we find the root?
if fc == 0:
print('f(c)==0')
break
elif np.sign(fa*fc) < 0:
b = c
fb = fc
else:
a = c
fa = fc
i += 1
xc = (a+b)/2.
return xc
In [6]:
## Finding a root of cos(x). What about if you change the interval?
#f = lambda x: np.cos(x)
## Another function
#f = lambda x: x**3-2*x**2+(4/3)*x-(8/27)
## Computing the cubic root of 7.
#f = lambda x: x**3-7
#bisect(f,0,2)
f = lambda x: x*np.exp(x)-3
#f2 = lambda x: np.cos(x)-x
bisect(f,0,3,tol=1e-13)
Out[6]:
It's very important to define a concept called convergence rate. This rate shows how fast the convergence of a method is at a specified point.
The convergence rate for the bisection is always 0.5 because this method uses the half of the interval for each iteration.
In [7]:
def cobweb(x,g=None):
min_x = np.amin(x)
max_x = np.amax(x)
plt.figure(figsize=(10,10))
ax = plt.axes()
plt.plot(np.array([min_x,max_x]),np.array([min_x,max_x]),'b-')
for i in np.arange(x.size-1):
delta_x = x[i+1]-x[i]
head_length = np.abs(delta_x)*0.04
arrow_length = delta_x-np.sign(delta_x)*head_length
ax.arrow(x[i], x[i], 0, arrow_length, head_width=1.5*head_length, head_length=head_length, fc='k', ec='k')
ax.arrow(x[i], x[i+1], arrow_length, 0, head_width=1.5*head_length, head_length=head_length, fc='k', ec='k')
if g!=None:
y = np.linspace(min_x,max_x,1000)
plt.plot(y,g(y),'r')
plt.title('Cobweb diagram')
plt.grid(True)
plt.show()
To learn about the Fixed-Point Iteration we will first learn about the concept of Fixed Point.
A Fixed Point of a function $g$ is a real number $r$, where $g(r) = r$
The Fixed-Point Iteration is based in the Fixed Point concept and works like this to find the root of a function:
\begin{equation} x_{0} = initial\_guess \\ x_{i+1} = g(x_{i})\end{equation}To find an equation's solution using this method you'll have to move around some things to rearrange the equation in the form $x = g(x)$. That way, you'll be iterating over the funcion $g(x)$, but you will not find $g$'s root, but $f(x) = g(x) - x$ (or $f(x) = x - g(x)$)'s root. In our following example, we'll find the solution to $f(x) = x - \cos(x)$ by iterating over the funcion $g(x) = \cos(x)$.
In [8]:
def fpi(g, x0, k, flag_cobweb=False):
x = np.empty(k+1)
x[0] = x0
error_i = np.nan
print(' i | x(i) | x(i+1) ||x(i+1)-x(i)| | e_i/e_{i-1}')
print('--------------------------------------------------------------')
for i in range(k):
x[i+1] = g(x[i])
error_iminus1 = error_i
error_i = abs(x[i+1]-x[i])
print('%2d | %.10f | %.10f | %.10f | %.10f' %
(i,x[i],x[i+1],error_i,error_i/error_iminus1))
if flag_cobweb:
cobweb(x,g)
return x[-1]
In [24]:
g = lambda x: np.cos(x)
fpi2(g, 2, 20, True)
Out[24]:
Let's quickly explain the Cobweb Diagram we have here. The blue line is the function $x$ and the red is the function $g(x)$. The point in which they meet is $g$'s fixed point. In this particular example, we start at $y = x = 1.5$ (the top right corner) and then we "jump" vertically to $y = \cos(1.5) \approx 0.07$. After this, we jump horizontally to $x = \cos(1.5) \approx 0.07$. Then, we jump again vertically to $y = \cos\left(\cos(1.5)\right) \approx 0.997$ and so on. See the pattern here? We're just iterating over $x = \cos(x)$, getting closer to the center of the diagram where the fixed point resides, in $x \approx 0.739$.
It's very important to mention that the algorithm will converge only if the rate of convergence $S < 1$, where $S = \left| g'(r) \right|$. If you want to use this method, you'll have to construct $g(x)$ starting from $f(x)$ accordingly. In this example, $g(x) = \cos(x) \Rightarrow g'(x) = -\sin(x)$ and $|-\sin(0.739)| \approx 0.67$.
In [22]:
g = lambda x: -(3/2)*x**2+(11/2)*x-2
gp = lambda x: -3*x+11/2
a=-1/2.7
g2 = lambda x: x+a*(x-g(x))
#x=np.linspace(2,3,100)
#plt.plot(x,gp(x),'-')
#plt.plot(x,gp(x)*0+1,'r-')
#plt.plot(x,gp(x)*0-1,'g-')
#plt.grid(True)
#plt.show()
fpi(g2, 2.45, 12, True)
Out[22]:
For this method, we want to iteratively find some function $f(x)$'s root, that is, the number $r$ for which $f(r) = 0$. The algorithm is as follows:
\begin{equation} x_0 = initial\_guess \end{equation}\begin{equation} x_{i+1} = x_i - \cfrac{f(x_i)}{f'(x_i)} \end{equation}which means that you won't be able to find $f$'s root if $f'(r) = 0$. In this case, you would have to use the modified version of this method, but for now let's focus on the unmodified version first. Newton's (unmodified) method is said to have quadratic convergence.
In [11]:
def newton_method(f, fp, x0, rel_error=1e-8, m=1):
#Initialization of hybrid error and absolute
hybrid_error = 100
error_i = np.inf
print('i | x(i) | x(i+1) | |x(i+1)-x(i)| | e_i/e_{i-1} | e_i/e_{i-1}^2')
print('----------------------------------------------------------------------------------------')
#Iteration counter
i = 1
while (hybrid_error > rel_error and hybrid_error < 1e12 and i < 1e4):
#Newton's iteration
x1 = x0-m*f(x0)/fp(x0)
#Checking if root was found
if f(x1) == 0.0:
hybrid_error = 0.0
break
#Computation of hybrid error
hybrid_error = abs(x1-x0)/np.max([abs(x1),1e-12])
#Computation of absolute error
error_iminus1 = error_i
error_i = abs(x1-x0)
#Increasing counter
i += 1
#Showing some info
print("%d | %.10f | %.10f | %.20f | %.10f | %.10f" %
(i, x0, x1, error_i, error_i/error_iminus1, error_i/(error_iminus1**2)))
#Updating solution
x0 = x1
#Checking if solution was obtained
if hybrid_error < rel_error:
return x1
elif i>=1e4:
print('Newton''s Method diverged. Too many iterations!!')
return None
else:
print('Newton''s Method diverged!')
return None
In [12]:
f = lambda x: np.sin(x)
fp = lambda x: np.cos(x) # the derivative of f
newton_method(f, fp, 3.1,rel_error=1e-14)
Out[12]:
In [32]:
f = lambda x: x**2
fp = lambda x: 2*x # the derivative of f
newton_method(f, fp, 3.1,rel_error=1e-2, m=2)
Out[32]:
https://en.wikipedia.org/wiki/Wilkinson%27s_polynomial
Final question: Why is the root far far away from $16$?
In [33]:
x = sym.symbols('x', reals=True)
W=1
for i in np.arange(1,21):
W*=(x-i)
W # Printing W nicely
Out[33]:
In [35]:
# Expanding the Wilkinson polynomial
We=sym.expand(W)
We
Out[35]:
In [36]:
# Just computiong the derivative
Wep=sym.diff(We,x)
Wep
Out[36]:
In [37]:
# Lamdifying the polynomial to be used with sympy
P=sym.lambdify(x,We)
Pp=sym.lambdify(x,Wep)
In [38]:
# Using scipy function to compute a root
root = optimize.newton(P,16)
print(root)
Build a FPI such that given $x$ computes $\displaystyle \frac{1}{x}$. Write down your solution below or go and see the solution
In [ ]:
print('Please try to think and solve before you see the solution!!!')
In [ ]:
g1 = lambda x: 1-x**3
g2 = lambda x: (1-x)**(1/3)
g3 = lambda x: (1+2*x**3)/(1+3*x**2)
fpi(g3, 0.5, 10, True)
In [ ]:
g1p = lambda x: -3*x**2
g2p = lambda x: -(1/3)*(1-x)**(-2/3)
g3p = lambda x: ((1+3*x**2)*(6*x**2)-(1+2*x**3)*6*x)/((1+3*x**2)**2)
r=0.6823278038280194
print(g3p(r))
In [18]:
def fpi2(g, x0, k, flag_cobweb=False):
x = np.empty(k+1)
x[0] = x0
error_i = np.inf
print(' i | x(i) | x(i+1) ||x(i+1)-x(i)| | e_i/e_{i-1} | e_i/e_{i-1}^2')
print('-----------------------------------------------------------------------------')
for i in range(k):
x[i+1] = g(x[i])
error_iminus1 = error_i
error_i = abs(x[i+1]-x[i])
print('%2d | %.10f | %.10f | %.10f | %.10f | %.10f' %
(i,x[i],x[i+1],error_i,error_i/error_iminus1,error_i/(error_iminus1**2)))
if flag_cobweb:
cobweb(x,g)
return x[-1]
Which function shows quadratic convergence? Why?
In [ ]:
g1 = lambda x: (4./5.)*x+1./x
g2 = lambda x: x/2.+5./(2*x)
g3 = lambda x: (x+5.)/(x+1)
fpi2(g1, 3.0, 30, True)
In [ ]:
# What is 'a'? Can we find another 'a'?
a = -3*(1.7**2)
print(a)
In [ ]:
f = lambda x: x**3-7
g = lambda x: f(x)/a+x
r=fpi(g, 1.7, 14, True)
print(f(r))
In [ ]:
f = lambda x: 8*x**4-12*x**3+6*x**2-x
fp = lambda x: 32*x**3-36*x**2+12*x-1
x = np.linspace(-1,1,1000)
plt.figure(figsize=(10,10))
plt.title('What are we seeing with the semiloigy plot? Is this function differentiable?')
plt.semilogy(x,np.abs(f(x)),'b-')
plt.semilogy(x,np.abs(fp(x)),'r-')
plt.grid()
plt.ylabel('$f(x)$',fontsize=16)
plt.xlabel('$x$',fontsize=16)
plt.show()
In [ ]:
r=newton_method(f, fp, 0.3, rel_error=1e-8, m=1)
print([r,f(r)])
# Is this showing quadratic convergence? If not, can you fix it?
In [20]:
# We are finding the 1/a
# Solution code:
a = 2.1
g = lambda x: 2*x-a*x**2
gp = lambda x: 2-2*a*x
r=fpi2(g, 0.7, 7, flag_cobweb=True)
print([r,1/a])
# Are we seeing quadratic convergence?
In [21]:
xx=np.linspace(0.2,0.8,1000)
plt.figure(figsize=(10,10))
plt.plot(xx,g(xx),'-',label=r'$g(x)$')
plt.plot(xx,gp(xx),'r-',label=r'$gp(x)$')
plt.plot(xx,xx,'g-',label=r'$x$')
plt.plot(xx,0*xx+1,'k--')
plt.plot(xx,0*xx-1,'k--')
plt.legend(loc='best')
plt.grid()
plt.show()
In [ ]:
g = lambda x: (x**2-1.)/2.
gh = lambda x: x+0.7*(x-g(x))
#fpi(g, 2, 15, True)
fpi2(gh, 0, 15, True)
In [ ]:
g1 = lambda x: np.cos(x)
g2 = lambda x: np.cos(x)**2
g3 = lambda x: np.sin(x)
g4 = lambda x: 1-5*x+(15/2)*x**2-(5/2)*x**3 # (0,1) y (1,2)
interact(lambda x0, n: fpi2(g4, x0, n, True),x0=widgets.FloatSlider(min=0, max=3, step=0.01, value=0), n=widgets.IntSlider(min=10, max=100, step=1, value=10))
In [ ]: