In [57]:
def f(x):
return x**3 - x +1
def df(x):
return 3*x**2 - 1
xmesh = np.linspace(-2, 2, 100)
pt.ylim([-3, 10])
pt.plot(xmesh, f(xmesh))
Out[57]:
import numpy as np
import scipy as sp
import matplotlib.pyplot as pt
import scipy.optimize as opt
def f(x):
return x**3 - x +1
def df(x):
return 3*x**2 - 1
xmesh = np.linspace(-2, 2, 100)
pt.ylim([-3, 10])
pt.plot(xmesh, f(xmesh))
iterates = [2, 1.5]
x_exact = opt.fsolve(f, -1)
print "Exact:", x_exact
e = list(np.abs(iterates-x_exact))
print "Errors:", e
# Evaluate this cell many times in-place (using Ctrl-Enter)
xlast, x = iterates[-2:] # grab last two iterates
slope = df(x) # Newton's method
#slope = ( (f(x)-f(xlast)) / (x-xlast) ) # Secant method
pt.plot(xmesh, f(xmesh))
pt.plot(xmesh, f(x) + slope*(xmesh-x))
pt.plot(x, f(x), "o")
pt.ylim([-3, 10])
pt.axhline(0, color="black")
# Compute approximate root
xnew = x - f(x) / slope
iterates.append(xnew)
print xnew
e.append(abs(xnew - x_exact))
for i, x in enumerate(iterates):
if i == 0:
continue
print "%3d %10g %10g | %10g %10g %10g"%(
i, x, e[i],
e[i]/e[i-1],
e[i]/e[i-1]**2,
e[i]/e[i-1]**1.62)