# Recursive Taylor Series Using Python's Sympy Library¶

Greetings everyone, I ran into some homework problems involving taking a multivariate Taylor expansion and noticed that Sympy's doesn't have this feature. They do have a series expansion at one point. However, things go sour for additional terms. I managed to write a multivariate Taylor series expanded about the origin up to the second order and provided the code below for those who are interested. If you see something that can be improved, let me know at ttucker@sdsu.edu. Thanks!

Getting started, The Taylor series for a function $f(x)$ centered at $x=a$

$$f(x) = f(a) + (x-a)\frac{\partial f(x)}{\partial x} + (x-a)\frac{\partial f(x)}{\partial x} + ... +$$ $$\frac{(x-a)^n}{n!}\frac{\partial^n f(x)}{\partial x^n}+ ... + \frac{(x-a)^{n+1}}{(n+1)!} \frac{\partial^{n+1} f(x)}{\partial x^{n+1}} + ...$$ $$= \sum_{j=0}^{\infty} \frac{(x-a)^j}{j!}\frac{\partial^j f(x)}{\partial x^j}$$

What about a multivariate function $v(\vec{x})$? Where $\vec{x} = [x_1, x_2, ..., x_n]$. Centered about $\vec{a}$ the second order Taylor expansion is $$g(\vec{x}) \approx g(\vec{a}) + \nabla g(\vec{a}) \cdot (\vec{x} - \vec{a}) + (\vec{x} - \vec{a})^T \cdot H(\vec{x}) \cdot (\vec{x} - \vec{a})$$ Where $\nabla g(a)$ is the gradient of g at point $\vec{a}$ and $H(\vec{x})$ is the Hessian matrix. Can we go higher? Of course we can, but this is where vector notation breaks down. Its a shame, because its so convinient to program in, but we can't expand past this term. The third way to write the Taylor expansion offers a way to extend the notation. $$g(a + \Delta\vec{x}) = \sum_j^\infty \frac{(\Delta\vec{x}\nabla)^jg(a)}{j!}$$ Where $$(\Delta\vec{x}\nabla)^j = \sum_{|a| \leq j}^n \frac{j!}{a!}\Delta\vec{x}^{\alpha}\partial^\alpha$$ For each variable $\alpha = x_1, ... x_n$. This is getting out of hand...however, this last equation offers us a way to write the Taylor series out recursively. I will show you with code below.

In [65]:
import sympy as sy
from sympy.matrices import *
x, t, a, h, k = sy.symbols('x t a h k')
x0, t0 = sy.symbols('x_0, t_0')

v = sy.Function('v')(x, t)
X = Matrix([x,t])
dX = Matrix([h,k])
sy.init_printing()


We are using sympy's matrix notation. Matrices are denoted with capital letters ($H$). I indicate Vectors by an arrow on top ($\vec{x}$). Sympy runs into trouble when you try substituting $x$ with $x-x_0$, so we are stuck with Taylor expanding about the origin.

In [66]:
dX

Out[66]:
$$\left[\begin{matrix}h\\k\end{matrix}\right]$$
In [67]:
X

Out[67]:
$$\left[\begin{matrix}x\\t\end{matrix}\right]$$

### The zeroth term¶

We kick off the Taylor series by adding the zeroth order term. Its just v calculated at the point $[0, 0]$.

In [68]:
TE = 0
v0 = v
for idx, x in enumerate(X):
v0 = v0.subs(x, 0)
TE += v0

In [69]:
TE

Out[69]:
$$v{\left (0,0 \right )}$$

### The first order terms¶

The first derivative portion is the gradient transpose multiplied by $\vec{\Delta x}$

$$\nabla v \cdot \vec{\Delta x}$$

In [70]:
gradV = Matrix([v.diff(xi,1) for xi in X])

In [71]:
gradV

Out[71]:
$$\left[\begin{matrix}\frac{\partial}{\partial x} v{\left (x,t \right )}\\\frac{\partial}{\partial t} v{\left (x,t \right )}\end{matrix}\right]$$
In [72]:
firstTerm = gradV.T*dX
TE += sum(firstTerm)

In [73]:
firstTerm[0]

Out[73]:
$$h \frac{\partial}{\partial x} v{\left (x,t \right )} + k \frac{\partial}{\partial t} v{\left (x,t \right )}$$

We will soon see that it's better to rewrite this as a function called applyDiff

In [74]:
def factorial(n):
if n <= 0:
return 1
else:
return n * factorial(n - 1)

In [75]:
tsOrderGroup = lambda f, X, n: dX.T * Matrix([f.diff(xi,1) for xi in X]) / factorial(n)
grad = tsOrderGroup(v, X, 1).T
fir = grad
sum(fir)

Out[75]:
$$h \frac{\partial}{\partial x} v{\left (x,t \right )} + k \frac{\partial}{\partial t} v{\left (x,t \right )}$$

### Second order terms¶

Next, we have to extend this to the Hessian matrix $H$. The second order term is

$$\frac{1}{2} \Delta \vec{x}^T H \Delta \vec{x}$$

In [76]:
sy.hessian(v, X)

Out[76]:
$$\left[\begin{matrix}\frac{\partial^{2}}{\partial x^{2}} v{\left (x,t \right )} & \frac{\partial^{2}}{\partial t\partial x} v{\left (x,t \right )}\\\frac{\partial^{2}}{\partial t\partial x} v{\left (x,t \right )} & \frac{\partial^{2}}{\partial t^{2}} v{\left (x,t \right )}\end{matrix}\right]$$
In [77]:
secondTerm = .5*dX.T*sy.hessian(v, X)*dX
TE += sum(secondTerm)
sy.simplify(sum(secondTerm))

Out[77]:
$$0.5 h^{2} \frac{\partial^{2}}{\partial x^{2}} v{\left (x,t \right )} + 1.0 h k \frac{\partial^{2}}{\partial t\partial x} v{\left (x,t \right )} + 0.5 k^{2} \frac{\partial^{2}}{\partial t^{2}} v{\left (x,t \right )}$$

Using our own function, we show the same thing.

In [78]:
gNext = Matrix([tsOrderGroup(dv, X, 2) for dv in grad * factorial(2 - 1)])
sy.simplify(sum(gNext))

Out[78]:
$$\frac{h^{2}}{2} \frac{\partial^{2}}{\partial x^{2}} v{\left (x,t \right )} + h k \frac{\partial^{2}}{\partial t\partial x} v{\left (x,t \right )} + \frac{k^{2}}{2} \frac{\partial^{2}}{\partial t^{2}} v{\left (x,t \right )}$$

### remaining terms with $\mathcal{O}(\cdot)$¶

We can continue expanding, or just denote the lowest order terms not showing in $\mathcal{O}(\cdot)$ notation.

In [79]:
gNext = Matrix([tsOrderGroup(dv, X, 3) for dv in gNext * factorial(3 - 1)])
TE += sum(gNext)
sy.simplify(sum(gNext))

Out[79]:
$$\frac{h^{3}}{6} \frac{\partial^{3}}{\partial x^{3}} v{\left (x,t \right )} + \frac{h^{2} k}{2} \frac{\partial^{3}}{\partial t\partial x^{2}} v{\left (x,t \right )} + \frac{h k^{2}}{2} \frac{\partial^{3}}{\partial t^{2}\partial x} v{\left (x,t \right )} + \frac{k^{3}}{6} \frac{\partial^{3}}{\partial t^{3}} v{\left (x,t \right )}$$

### remaining terms with $\mathcal{O}(\cdot)$¶

We can continue expanding, or just denote the lowest order terms not showing in $\mathcal{O}(\cdot)$ notation.

In [80]:
TE += sy.O(h**4) + sy.O(k**4)

In [81]:
sy.simplify(TE)

Out[81]:
$$1.0 v{\left (0,0 \right )} + 1.0 k \frac{\partial}{\partial t} v{\left (x,t \right )} + 0.5 k^{2} \frac{\partial^{2}}{\partial t^{2}} v{\left (x,t \right )} + 0.166666666666667 k^{3} \frac{\partial^{3}}{\partial t^{3}} v{\left (x,t \right )} + 1.0 h \frac{\partial}{\partial x} v{\left (x,t \right )} + 1.0 h k \frac{\partial^{2}}{\partial t\partial x} v{\left (x,t \right )} + 0.5 h k^{2} \frac{\partial^{3}}{\partial t^{2}\partial x} v{\left (x,t \right )} + 0.5 h^{2} \frac{\partial^{2}}{\partial x^{2}} v{\left (x,t \right )} + 0.5 h^{2} k \frac{\partial^{3}}{\partial t\partial x^{2}} v{\left (x,t \right )} + 0.166666666666667 h^{3} \frac{\partial^{3}}{\partial x^{3}} v{\left (x,t \right )} + \mathcal{O}\left(k^{4}\right) + \mathcal{O}\left(h^{4}\right)$$

We can remove the $\mathcal{O}(\cdot)$ any time with the removeO() method.

In [82]:
TE.removeO()

Out[82]:
$$\frac{h}{6} \left(h \left(h \frac{\partial^{3}}{\partial x^{3}} v{\left (x,t \right )} + k \frac{\partial^{3}}{\partial t\partial x^{2}} v{\left (x,t \right )}\right) + k \left(h \frac{\partial^{3}}{\partial t\partial x^{2}} v{\left (x,t \right )} + k \frac{\partial^{3}}{\partial t^{2}\partial x} v{\left (x,t \right )}\right)\right) + h \left(0.5 h \frac{\partial^{2}}{\partial x^{2}} v{\left (x,t \right )} + 0.5 k \frac{\partial^{2}}{\partial t\partial x} v{\left (x,t \right )}\right) + h \frac{\partial}{\partial x} v{\left (x,t \right )} + \frac{k}{6} \left(h \left(h \frac{\partial^{3}}{\partial t\partial x^{2}} v{\left (x,t \right )} + k \frac{\partial^{3}}{\partial t^{2}\partial x} v{\left (x,t \right )}\right) + k \left(h \frac{\partial^{3}}{\partial t^{2}\partial x} v{\left (x,t \right )} + k \frac{\partial^{3}}{\partial t^{3}} v{\left (x,t \right )}\right)\right) + k \left(0.5 h \frac{\partial^{2}}{\partial t\partial x} v{\left (x,t \right )} + 0.5 k \frac{\partial^{2}}{\partial t^{2}} v{\left (x,t \right )}\right) + k \frac{\partial}{\partial t} v{\left (x,t \right )} + v{\left (0,0 \right )}$$

I will now wrap this all up into a general Taylor series function.

In [83]:
def myTaylorSeries(v, nOrder, X, dX):
TE = 0
v0 = v
for idx, x in enumerate(X):
v0 = v0.subs(x, 0)
TE += sy.O(x**(nOrder+1))
TE += v0
tsOrderGroup = lambda f, X, n: dX.T * Matrix([f.diff(xi,1) for xi in X])
grad = tsOrderGroup(v, X, 1).T
TE += sum(grad)
gNext = grad
for nord in range(2, nOrder+1):
gNext = Matrix([tsOrderGroup(dv, X, nord) for dv in gNext])
TE += sum(gNext)/factorial(nord)
return TE

In [84]:
sy.simplify(myTaylorSeries(v, 4, X, dX))

Out[84]:
$$v{\left (0,0 \right )} + k \frac{\partial}{\partial t} v{\left (x,t \right )} + \frac{k^{2}}{2} \frac{\partial^{2}}{\partial t^{2}} v{\left (x,t \right )} + \frac{k^{3}}{6} \frac{\partial^{3}}{\partial t^{3}} v{\left (x,t \right )} + \frac{k^{4}}{24} \frac{\partial^{4}}{\partial t^{4}} v{\left (x,t \right )} + h \frac{\partial}{\partial x} v{\left (x,t \right )} + h k \frac{\partial^{2}}{\partial t\partial x} v{\left (x,t \right )} + \frac{h k^{2}}{2} \frac{\partial^{3}}{\partial t^{2}\partial x} v{\left (x,t \right )} + \frac{h k^{3}}{6} \frac{\partial^{4}}{\partial t^{3}\partial x} v{\left (x,t \right )} + \frac{h^{2}}{2} \frac{\partial^{2}}{\partial x^{2}} v{\left (x,t \right )} + \frac{h^{2} k}{2} \frac{\partial^{3}}{\partial t\partial x^{2}} v{\left (x,t \right )} + \frac{h^{2} k^{2}}{4} \frac{\partial^{4}}{\partial t^{2}\partial x^{2}} v{\left (x,t \right )} + \frac{h^{3}}{6} \frac{\partial^{3}}{\partial x^{3}} v{\left (x,t \right )} + \frac{h^{3} k}{6} \frac{\partial^{4}}{\partial t\partial x^{3}} v{\left (x,t \right )} + \frac{h^{4}}{24} \frac{\partial^{4}}{\partial x^{4}} v{\left (x,t \right )} + \mathcal{O}\left(x^{5}\right) + \mathcal{O}\left(t^{5}\right)$$

Let's try a 4d system.

In [85]:
x, y, z, t, h, l, m, k = sy.symbols('x, y, z, t, h, l, m, k')
S = Matrix([x, y, z, t])
dS = Matrix([h, l, m, k])
f = sy.Function('f')(x, y, z, t)
fTS = myTaylorSeries(f, 2, S, dS)

In [86]:
sy.simplify(fTS)

Out[86]:
$$f{\left (0,0,0,0 \right )} + m \frac{\partial}{\partial z} f{\left (x,y,z,t \right )} + \frac{m^{2}}{2} \frac{\partial^{2}}{\partial z^{2}} f{\left (x,y,z,t \right )} + l \frac{\partial}{\partial y} f{\left (x,y,z,t \right )} + l m \frac{\partial^{2}}{\partial y\partial z} f{\left (x,y,z,t \right )} + \frac{l^{2}}{2} \frac{\partial^{2}}{\partial y^{2}} f{\left (x,y,z,t \right )} + k \frac{\partial}{\partial t} f{\left (x,y,z,t \right )} + k m \frac{\partial^{2}}{\partial t\partial z} f{\left (x,y,z,t \right )} + k l \frac{\partial^{2}}{\partial t\partial y} f{\left (x,y,z,t \right )} + \frac{k^{2}}{2} \frac{\partial^{2}}{\partial t^{2}} f{\left (x,y,z,t \right )} + h \frac{\partial}{\partial x} f{\left (x,y,z,t \right )} + h m \frac{\partial^{2}}{\partial x\partial z} f{\left (x,y,z,t \right )} + h l \frac{\partial^{2}}{\partial x\partial y} f{\left (x,y,z,t \right )} + h k \frac{\partial^{2}}{\partial t\partial x} f{\left (x,y,z,t \right )} + \frac{h^{2}}{2} \frac{\partial^{2}}{\partial x^{2}} f{\left (x,y,z,t \right )} + \mathcal{O}\left(z^{3}\right) + \mathcal{O}\left(y^{3}\right) + \mathcal{O}\left(x^{3}\right) + \mathcal{O}\left(t^{3}\right)$$

# Summary:¶

I'm thrilled with this function. Taylor series crop up everywhere. They are used to derive new numerical schemes; they are used for analysis, they are used in derivations. And they are a pain in the neck to expand in higher dimensions until now. Please feel free to use this to your own devices. Thanks for reading!