In [1]:

%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
import sympy as sy
import cmath as cm
from scipy.misc import derivative

Newton’s Method¶

Newton’s method is a way to use the tangent line to approximate zeros of functions. A simple example to motivate the idea is to consider the example of the function \(f(x) = x^2 - 1\). Obviously this function has zeros at \(x = -1\) and \(x = 1\).

Below, we plot the function, and a line tangent to the graph of the function at the points \((1.1, f(1.1))\) and \((0.9, f(0.9))\).

In [2]:

def f(x):
    return (x**2 - 1)

def tan_line(a):
    return derivative(f, a)*(x - a) + f(a)

In [3]:

x = np.linspace(-2,2,100)

In [4]:

plt.plot(x, f(x))
plt.plot(x, tan_line(1.1))
plt.axhline(color = 'black')
plt.axvline(color = 'black')
plt.title("Tangent at x=1.1")

Out[4]:

<matplotlib.text.Text at 0x10f289978>

In [11]:

from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes
from mpl_toolkits.axes_grid1.inset_locator import mark_inset

fig, ax = plt.subplots(figsize = (10, 7)) # create a new figure with a default 111 subplot
ax.plot(x, f(x))
ax.plot(x, tan_line(a = 1.2))
plt.ylim(-2,2)
ax.axhline(color = 'black')
axins = zoomed_inset_axes(ax, 3.5, loc=2)
axins.plot(x, f(x))
axins.plot(x, tan_line(a = 1.2))


x1, x2, y1, y2 = 0.8,1.2, -0.2, 0.2
axins.set_xlim(x1, x2)
axins.set_ylim(y1, y2)
plt.yticks(visible=False)
plt.xticks(visible=False)
mark_inset(ax, axins, loc1=1, loc2=3, fc="none", ec="0.5")
axins.axhline(color = 'black')

Out[11]:

<matplotlib.lines.Line2D at 0x1127fd160>

In [5]:

def df(x):
    dx = 0.000001
    return (f(x + dx) - f(x))/dx

In [12]:

def N(x):
    return x - f(x)/df(x)

In [13]:

N(1)

Out[13]:

1.0

In [14]:

l = [1.3]
for i in range(8):
    z = N(l[i])
    l.append(z)

In [15]:

Out[15]:

[1.3,
0346154866690602,
0005790875782625,
0000001678633645,
000000000000098,
0,
0,
0,
0]

In [16]:

N(1.3478260274025091)

Out[16]:

1.0448808843824686

In [17]:

def f(x):
    return x**3 - x

In [18]:

l = [1+.5j]

In [19]:

for i in range(120000):
    z = N(l[i])
    l.append(z)

In [20]:

X = [x.real for x in l]
Y = [x.imag for x in l]

In [21]:

Out[21]:

[1.0,
8402370247174764,
9215077851434575,
9920869754787999,
9992702709978157,
0000004025142377,
9999999999989508,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
 ...]

In [22]:

def f(x):
    return x**2 - 1

from scipy.misc import derivative

x = np.linspace(-2, 2, 1000)
df = derivative(f, x)

In [24]:

plt.figure()
plt.plot(df)
plt.plot(f(x))

Out[24]:

[<matplotlib.lines.Line2D at 0x1075716d8>]

In [28]:

def N(x):
    return x - f(x)/derivative(f, x)

In [29]:

N(1)

Out[29]:

1.0

In [30]:

N(0)

/Users/NYCMath/anaconda/lib/python3.5/site-packages/ipykernel/__main__.py:2: RuntimeWarning: divide by zero encountered in double_scalars
  from ipykernel import kernelapp as app

Out[30]:

inf

In [31]:

N(.4)

Out[31]:

1.4500000000000002

In [34]:

plt.figure()
plt.plot(x, f(x))
plt.axhline(color = 'black')
plt.axvline(color = 'black')

Out[34]:

<matplotlib.lines.Line2D at 0x11397cf28>

In [33]:

N(-1)

Out[33]:

-1.0

In [35]:

N(N(-1))

Out[35]:

-1.0

In [36]:

N(-1.2)

Out[36]:

-1.0166666666666666

In [37]:

N(N(-1.2))

Out[37]:

-1.000136612021858

In [38]:

def ns(x0, N):
    for i in range(N):
        start = [x0]

In [39]:

def f(x):
    return (x-2)**2 - 2

x = np.linspace(0, 5, 100)

In [42]:

plt.figure()
plt.plot(x, f(x))
plt.axhline(color = 'black')
plt.axvline(color = 'black')

Out[42]:

<matplotlib.lines.Line2D at 0x113c6c2b0>

In [43]:

def df(a):
    return derivative(f, a)*(x - a) + f(a)

In [49]:

fig, ax = plt.subplots()
import matplotlib.patches as patches
plt.plot(x, f(x))
plt.plot(1.4, f(1.4), 'o', markersize = 10, alpha = 0.6)
plt.plot(x, df(a = 1.3))
plt.plot(x, df(a = 1.2))
plt.plot(x, df(a = 1.1))
plt.axhline(color = 'black')
plt.axvline(color = 'black')
plt.xlim(0,1.5)
plt.ylim(f(1.4)-0.2, 2)


plt.plot(x, df(a = 1.4))

Out[49]:

[<matplotlib.lines.Line2D at 0x114dd3828>]

In [50]:

def f(x):
    return x**2

x = np.linspace(-2,2,100)
width = (2-(-2))/100
difs = np.diff(f(x))
plt.figure(figsize=(10, 7))
plt.subplot(2, 2, 1)
plt.plot(x, f(x), '-o', markersize = 1)
plt.title("Function")
plt.subplot(2, 2, 2)
plt.plot(x[1:], difs, '-o', markersize = 1)
plt.axhline(color = 'black')
plt.title("Differences")

plt.subplot(2, 2, 3)
plt.plot(x,f(x))
plt.title("Function")
plt.subplot(2,2,4)
plt.plot(x, derivative(f, x))
plt.axhline(color = 'black')
plt.title("Derivative")

Out[50]:

<matplotlib.text.Text at 0x1150f9cf8>

In [ ]:

plt.plot(f(x), animated = True)

In [51]:

import numpy as np
import matplotlib.pyplot as plt

from matplotlib import animation, rc
from IPython.display import HTML

In [52]:

# First set up the figure, the axis, and the plot element we want to animate

fig, ax = plt.subplots()

ax.set_xlim(( 0, 2))
ax.set_ylim((-2, 2))

line, = ax.plot([], [], lw=2)

In [53]:

def init():
    line.set_data([], [])
    return (line,)

In [54]:

def animate(i):
    x = np.linspace(0, 2, 1000)
    y = np.sin(2 * np.pi * (x - 0.01 * i))
    line.set_data(x, y)
    return (line,)

In [55]:

anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=100, interval=20, blit=True)

HTML(anim.to_html5_video())

Out[55]:

In [56]:

rc('animation', html='html5')
anim

Out[56]: