import numpy as np
import matplotlib.pyplot as plt
%matplotlib notebook
def f(x):
return x**3 - 4*x**2 - 5.
xx = np.linspace(-1, 5, 2000)
plt.figure(figsize=(4,4))
plt.plot(xx, f(xx));
plt.grid('on');
plt.axhline(0, color='k', lw=0.75);
f(4.273748)
-1.4141239702780695e-05
def phi(x):
return np.cos(x)
x = np.linspace(0,1.5,200);
plt.figure(figsize=(4,4))
plt.plot(x, phi(x)); plt.axis('equal'); plt.grid('on');
plt.plot(x,x,'r');
def FixedPointIteration(x0, phi):
x_history = []
phi_history = []
for idx in range(10):
x_history.append(x0)
phi_history.append(phi(x0))
x0 = phi(x0)
x_history.append(x0)
phi_history.append(x0)
print(x0)
return x_history, phi_history
xh, phih = FixedPointIteration(0.25, phi)
0.9689124217106447 0.5661963244509127 0.8439474651598078 0.6645181803460104 0.7872140049933128 0.7058216222217144 0.7610788557592939 0.7240923421872705 0.7491010136260781 0.7323013570816245
plt.figure(figsize=(8,8))
plt.plot(x, phi(x)); plt.axis('equal'); plt.grid('on');
plt.plot(x,x,'r');
plt.plot(xh, phih,'g')
plt.plot(xh[0],phih[0],'go')
[<matplotlib.lines.Line2D at 0x1205a8490>]
def phi(x):
return x - (x**2 - 2.)
#%pylab auto
x = np.linspace(0,2,200);
plt.figure(figsize=(4,4))
plt.plot(x, phi(x));
plt.axis('equal'); plt.grid('on');
plt.plot(x,x,'r');
#plt.axis([x[0],x[-1],x[0],x[-1]])
plt.show()
xh, phih = FixedPointIteration(0.25, phi)
2.1875 -0.59765625 1.0451507568359375 1.9528106523212045 0.13934120850203602 2.119925236115228 -0.3741577706029773 1.4858481920944326 1.2781033421441386 1.6445551889441217
plt.figure(figsize=(4,4))
plt.plot(x, phi(x));
plt.axis('equal'); plt.grid('on');
plt.plot(x,x,'r');
#plt.axis([x[0],x[-1],x[0],x[-1]])
plt.plot(xh, phih,'g')
plt.plot(xh[0],phih[0],'go')
plt.show()
# Let's illustrate what Newton's method does.
def f(x):
return x**2 - 2.
def fp(x):
return 2.*x
x = np.linspace(0,2,200)
plt.figure(figsize=(8,8))
plt.plot(x,np.zeros(x.shape),'k', x, f(x)); plt.grid('on');
a = plt.axis()
x0 = 1.; f0 = f(x0);
m = fp(x0)
f_left = f0 + m*(a[0]-x0); f_right = f0 + m*(a[1]-x0);
plt.plot(a[:2], [f_left, f_right], 'r');
plt.axis(a);
plt.plot([x0], [f0], 'ro');
x1 = x0 - f0/m
plt.plot([x1], [0], 'ro');
# Now let's try it on phi(x) = cos(x)
def phi_Newton(x):
return x - (-np.cos(x))/np.sin(x)
xh, phih = FixedPointIteration(0.8, phi_Newton)
1.7712146006504743 1.5680690668947184 1.5707963335566546 1.5707963267948966 1.5707963267948966 1.5707963267948966 1.5707963267948966 1.5707963267948966 1.5707963267948966 1.5707963267948966
plt.figure()
plt.plot(x,x,'r', x, x+1./np.tan(x))
plt.axis([0.5, 2, 0.5, 2])
plt.plot(xh, phih,'g')
plt.plot(xh[0],phih[0],'go');
/var/folders/97/ppgk101x3214hn7bzlb85jk00000gn/T/ipykernel_54967/2142578811.py:2: RuntimeWarning: divide by zero encountered in divide plt.plot(x,x,'r', x, x+1./np.tan(x))
kk = np.arange(1,10)
c = 0.8
d = 0.5
linear = []
for k in kk:
linear.append(c*d**1)
d = linear[-1]
d = 0.5
quadratic = []
for k in kk:
quadratic.append(c*d**2)
d = quadratic[-1]
plt.figure(figsize=(4,4))
plt.plot(kk, linear, 'o');
plt.plot(kk, quadratic, 'o');
plt.ylabel(r'$|x_{k+1} - x_k|$'); plt.xlabel(r'$k$');
plt.legend(['Linear', 'Quadratic']);
plt.tight_layout();