Tuesday, February 8, 2011

Plotting the Taylor Series for sine

One day about a year ago I was in a chain bookstore killing time, not really expecting to find much, and I came across a Calculus textbook which has become one of my favorites. The author is Morris Kline, and the title, Calculus, An intuitive and Physical Approach, tells most of the story. It's a minor miracle that the book was in the store---it had obviously been there for years, but was too big for them to discard, or something.

Anyway, I was browsing Kline's section on Taylor Series last night and it made me think of a little project. Before I get to that, if you don't know, the Taylor Series for sin(x), cos(x) and ex can be used to derive Euler's famous formula, which we touched on briefly a while back (here and here).

The project: there is a terrific graphic in the wikipedia article showing what a good approximation one can get with only a few terms from the series. So that's our goal for today, to reproduce the plot using matplotlib.

The Taylor series for f(x) when x is near a:

f(a) + f1(a)/1!(x-a) + f2(a)/2!(x-a)2 + .. 

where fn is the nth derivative of f

If a = 0 this is the Maclaurin series

f(0) + f1(0)/1!(x) + f2(0)/2!(x)2 + ..

For sin(0) and cos(0) the sequential derivatives form a repeating series (with period equal to 4) and every other term equal to 0. If f(x) is sin(x):

sin(0) + cos(0)/1!(x) - sin(0)/2!(x)2 - cos(0)/3!(x)3
+ sin(0)4!(x)4 + ..

All the sin(0) terms are equal to 0, and cos(0) is equal to 1:

  = x - x3/3! + x5/5! - x7/7! ..

To make the plot, we define a "factory" function that returns functions for the various degrees:

def func(n):
def f(x):
v = x**n / factorial(n)
if n % 4 == 3:
v *= -1
return v
return f

which we will call with odd values of n. Then the only trick to the plot is to accumulate the sum. We plot the sine curve last so that it lies on top of the others. Sure enough, with 8 terms we can cover roughly -2π < x < 2π.

from __future__ import division
from math import pi, factorial
import matplotlib.pyplot as plt
import numpy as np

X = np.arange(-3*pi, 3*pi, 0.001)

def func(n):
def f(x):
v = x**n / factorial(n)
if n % 4 == 3:
v *= -1
return v
return f

cL = ['r','g','b','purple','cyan',
D = dict()

for i in range(1,16,2):
c = cL.pop(0)
D[i] = func(i)(X)
j = i-2
if j > 0:
D[i] += D[j]


ax = plt.axes()