## Taylor expansions¶

We first import some prerequisite modules, most importantly Sympy and Matplotlib.

In [1]:
def uniq(seq):
# order preserving
seen = {}
result = []
for item in seq:
if item in seen:
continue
seen[item] = 1
result.append(item)
return result

In [2]:
import sympy as sm
from sympy.abc import x
from matplotlib.colors import cnames

In [3]:
import matplotlib
from matplotlib.colors import rgb2hex

from numpy import arange

In [4]:
from sympy import Symbol,plot
import matplotlib.pyplot as plt

def move_sympyplot_to_axes(p, ax):
backend = p.backend(p)
backend.ax = ax
backend.process_series()
backend.ax.spines['right'].set_color('none')
backend.ax.spines['bottom'].set_position('zero')
backend.ax.spines['top'].set_color('none')
plt.close(backend.fig)


The following routine will plot the function f together with its Taylor polynomials up to the degree degree, in the square box with sides [-A,A]. We will use it several times for different functions f.

In [5]:
def taylorplot(f, degree, A):
L = [f.series(n=i).removeO() for i in range (1, degree+2)] # +1 for the offset of range, another +1 for the offset of .series
L = uniq(L)
L.insert(0,f)
# Pick colors from the Blues colorscheme
a = arange(0,1,1/len(L))
c = matplotlib.cm.Blues(a)[:,:3]

# Plot the original function in red, then add the graphs of Taylor polynomials in different shades of blue
p1 = sm.plot(L[0], xlim=[-A,A], ylim=[-A,A], line_color = 'r',  show=False)
for i in range(1,len(L)):
p1.extend(sm.plot(L[i], xlim=[-A,A], ylim=[-A,A], line_color = rgb2hex(c[i,:]), show=False))
p1.legend = False
# p1.show()

# Plot the original function alone
p2 = sm.plot(L[0], xlim=[-A,A], ylim=[-A,A], line_color = 'r', show=False)
p2.legend = False
p2.grid = True
# p2.show()

# Generate the legend by LaTeXing the expressions for the functions involved.
l =[]
for i in range(len(L)):
l.append(sm.latex(L[i], mode='inline'))

# %matplotlib gtk3 # uncomment to use the external plot window
# Combine the two resulting plots, add legend
fig, (ax1,ax2) = plt.subplots(ncols=2,figsize=(20,10))
move_sympyplot_to_axes(p1, ax1)
move_sympyplot_to_axes(p2, ax2)
fig.legend(l,loc='lower center', bbox_to_anchor=(0.5, 0.0),
plt.show()


Let us first plot the partial sums of the Maclaurin series of the Bessel function $J_0$, up to the degree 14. We will restrict the graphs to the box [-10,10] by [-10,10]:

In [6]:
B = sm.besselj(0,x)
taylorplot(B, 14, 10)


In these graphs we show the Taylor polynomials (centered at $a=0$ here), using the darker shade to denote to the higher degree polynomial $T_n(x)$. The plot on the right shows $J_0(x)$ itself.

Next, let us do the same with $\sin(x)$. We plot its Taylor polynomials at 0 up to degree 13 in the box [-10,10] by [-10,10].

In [7]:
S = sm.sin(x)
taylorplot(S, 13,10)


Finally, let us look at two examples in which the Maclaurin series of a function does not converge for all x. The first example is $1/(1-x)$; we will plot its Taylor polynomials up to degree 10 in the square box [-6,6] by [-6,6].

In [8]:
f1 = 1/(1-x)
taylorplot(f1,10,6)


Observe how the polynomials $T_n$ diverge away from the graph of the function f1 when $x=-1$, even though the function does not have a singularity there. This is of course because the radius of convergence of the geometric series $$\sum_{n=0}^\infty x^n$$ is $1$. Since the interval of convergence of a series must be symmetric, divergence of the series at $x = 1$, where the function f1 has a singularity, means that the series must also diverge at $x = -1$.

The following example is even more subtle. The function we look at here, f2$= 1/(1+x^2)$, does not exhibit any singularities at all! On the other hand, its Maclaurin series $$\sum_{n=0}^\infty\, (-1)^n x^{2n}$$ diverges as soon as $|x|\geq 1$. Let us plot the Taylor polynomials at $a = 0$ up to the degree $n = 20$.

In [9]:
f2 = 1/(1+x**2)
taylorplot(f2,20,6)


Outside of Calc2: It turns out, to understand the reason for such behavior of this series, one has to extend the function f2$=1/(1+x^2)$ to complex numbers, where for $x = i$, our function f2 has a singularity (the complex number $i$ is such that $i^2+1 = 0$, turning the denominator of f2 into $0$). In the complex plane, series converge on discs instead of intervals, so because the function f2 has a singularity at the point $0+i\,1$, its series centered at $a=0$ will have the radius of convergence equal to $1$. We show this relation in the following graph. In it, a point $(x,y)$ in the plane stands for the complex number $z = x+ iy$.

In [10]:
    fig, ax = plt.subplots(figsize=(10,10))
circle =  matplotlib.patches.Circle ((0, 0), 1, edgecolor='none',alpha=0.5)
p = matplotlib.collections.PatchCollection([circle], alpha=0.4)
ax.scatter([0.0],[1.0],marker='o',c='r',s=40)
ax.annotate('Singularity of f2 at $0+i$', xy=(0.1, 1.05), xytext=(1, 1.5),
arrowprops=dict(arrowstyle="->")
)
ax.annotate('For complex $z= x + i\,y$\nthat lie inside this disc,\nthe series converges.', xy=(0.5, -0.5), xytext=(0.5, -1.5),
arrowprops=dict(arrowstyle="->")
)
ax.annotate('The series will diverge at the\n points of the disc boundary', xy=(-1.02, -0.13), xytext=(-2, -1),
arrowprops=dict(arrowstyle="->")
)
# This adjusts the position of the axes
ax.spines['left'].set_position('center')
ax.spines['bottom'].set_position('center')
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')

plt.ylim(-2,2)
plt.xlim(-2,2)
plt.show()