Sympy for computing Taylor series¶

(first part from Day 5)

In [10]:
import sympy as sp
In [11]:
t = sp.symbols('t')
h = sp.symbols('h')
In [12]:
(t+h)**2
Out[12]:
$\displaystyle \left(h + t\right)^{2}$
In [13]:
sp.expand((t+h)**2)
Out[13]:
$\displaystyle h^{2} + 2 h t + t^{2}$
In [14]:
sp.diff((t+h)**2, t)
Out[14]:
$\displaystyle 2 h + 2 t$
In [15]:
t,h = sp.symbols('t,h')
In [16]:
sp.expand((t+h)**3)
Out[16]:
$\displaystyle h^{3} + 3 h^{2} t + 3 h t^{2} + t^{3}$
In [17]:
alpha,Omega = sp.symbols('alpha,Omega')
alpha
Out[17]:
$\displaystyle \alpha$
In [18]:
Omega
Out[18]:
$\displaystyle \Omega$

We also need function objects¶

In [19]:
u = sp.Function('u')
u(t)
Out[19]:
$\displaystyle u{\left(t \right)}$
In [20]:
sp.diff(u(t),t)
Out[20]:
$\displaystyle \frac{d}{d t} u{\left(t \right)}$
In [21]:
sp.diff(u(alpha,t),t,t)
Out[21]:
$\displaystyle \frac{\partial^{2}}{\partial t^{2}} u{\left(\alpha,t \right)}$

Let's try to solve u' = pu, v'=-qv¶

where p and q depend on t

In [22]:
import sympy as sp
t,h =sp.symbols('t,h')
u = sp.Function('u')
v = sp.Function('v')
p = sp.Function('p')
q = sp.Function('q')
In [23]:
uprime = sp.diff(u(t),t)  # shorthand for the first derivative
vprime = sp.diff(v(t),t)

de = {uprime:p(t)*v(t),  vprime:-q(t)*u(t)}  # allow replacement according to the DE
de
Out[23]:
{Derivative(u(t), t): p(t)*v(t), Derivative(v(t), t): -q(t)*u(t)}

Continuing on Day 6:

In [ ]:
In [24]:
# python dict: an associative array linking "keys" (indices) to "values"
d = {'apple':'round fruit','zebra':'striped animal'}
d
Out[24]:
{'apple': 'round fruit', 'zebra': 'striped animal'}
In [25]:
d['zebra']
Out[25]:
'striped animal'
In [26]:
# sympy substitution

expr = t**2+h
expr
Out[26]:
$\displaystyle h + t^{2}$
In [27]:
expr.subs( {t:7}  )
Out[27]:
$\displaystyle h + 49$
In [28]:
expr.subs( {t:7,h:10}  )
Out[28]:
$\displaystyle 59$
In [29]:
mylist = [5,6,7]
mylist[0]
Out[29]:
5
In [30]:
uders = [u,de[uprime]]   # start with this, then append higher derivatives
vders = [v,de[vprime]]
In [31]:
uders
Out[31]:
[u, p(t)*v(t)]
In [32]:
# next element in uders will be the second derivative u''
sp.diff(uders[-1])
Out[32]:
$\displaystyle p{\left(t \right)} \frac{d}{d t} v{\left(t \right)} + v{\left(t \right)} \frac{d}{d t} p{\left(t \right)}$
In [33]:
# we also want to substitute for v' using the de
sp.diff(uders[-1]).subs(de)
Out[33]:
$\displaystyle - p{\left(t \right)} q{\left(t \right)} u{\left(t \right)} + v{\left(t \right)} \frac{d}{d t} p{\left(t \right)}$
In [34]:
uders = [u,de[uprime]]   # start with this, then append higher derivatives
vders = [v,de[vprime]]
for k in range(3):
    uders.append(  sp.diff(uders[-1]).subs(de)   )
    vders.append(  sp.diff(vders[-1]).subs(de)   )
In [35]:
uders
Out[35]:
[u,
 p(t)*v(t),
 -p(t)*q(t)*u(t) + v(t)*Derivative(p(t), t),
 -p(t)**2*q(t)*v(t) - p(t)*u(t)*Derivative(q(t), t) - 2*q(t)*u(t)*Derivative(p(t), t) + v(t)*Derivative(p(t), (t, 2)),
 p(t)**2*q(t)**2*u(t) - 2*p(t)**2*v(t)*Derivative(q(t), t) - 4*p(t)*q(t)*v(t)*Derivative(p(t), t) - p(t)*u(t)*Derivative(q(t), (t, 2)) - 3*q(t)*u(t)*Derivative(p(t), (t, 2)) - 3*u(t)*Derivative(p(t), t)*Derivative(q(t), t) + v(t)*Derivative(p(t), (t, 3))]

Now let's choose a specific DE, i.e. choose specific functions p,q¶

In [36]:
de = {uprime:sp.sin(t)*v(t),  vprime:-t*u(t)}  # allow replacement according to the DE
de
Out[36]:
{Derivative(u(t), t): v(t)*sin(t), Derivative(v(t), t): -t*u(t)}
In [37]:
uders = [u(t),de[uprime]]   # start with this, then append higher derivatives
vders = [v(t),de[vprime]]

for k in range(3): # 3 is how many higher-order derivatives I want to append (could be as large a number as you like)
    uders.append(  sp.diff(uders[-1]).subs(de)   )
    vders.append(  sp.diff(vders[-1]).subs(de)   )
In [38]:
uders  # our list of [u,u',u'',u''', ...]
Out[38]:
[u(t),
 v(t)*sin(t),
 -t*u(t)*sin(t) + v(t)*cos(t),
 -2*t*u(t)*cos(t) - t*v(t)*sin(t)**2 - u(t)*sin(t) - v(t)*sin(t),
 t**2*u(t)*sin(t)**2 + 3*t*u(t)*sin(t) - 4*t*v(t)*sin(t)*cos(t) - 3*u(t)*cos(t) - 2*v(t)*sin(t)**2 - v(t)*cos(t)]
In [39]:
# python sum
sum([1,2,3])
Out[39]:
6
In [41]:
# "list comprehension" for building a list whose elements are parametrized by another list-like thing
[i**2 for i in range(4)]
Out[41]:
[0, 1, 4, 9]
In [42]:
# now build the Taylor polynomial for u and for v
# Here are the first 5 terms:
[ h**i/sp.factorial(i)*uders[i]  for i in range(len(uders))]   # just need to sum these to get the Taylor polynomial
Out[42]:
[u(t),
 h*v(t)*sin(t),
 h**2*(-t*u(t)*sin(t) + v(t)*cos(t))/2,
 h**3*(-2*t*u(t)*cos(t) - t*v(t)*sin(t)**2 - u(t)*sin(t) - v(t)*sin(t))/6,
 h**4*(t**2*u(t)*sin(t)**2 + 3*t*u(t)*sin(t) - 4*t*v(t)*sin(t)*cos(t) - 3*u(t)*cos(t) - 2*v(t)*sin(t)**2 - v(t)*cos(t))/24]
In [43]:
# Here are the Taylor polynomials for u(t+h), v(t+h):
u_poly = sum([ h**i/sp.factorial(i)*uders[i]  for i in range(len(uders))]  )
v_poly = sum([ h**i/sp.factorial(i)*vders[i]  for i in range(len(vders))]  )
u_poly
Out[43]:
$\displaystyle \frac{h^{4} \left(t^{2} u{\left(t \right)} \sin^{2}{\left(t \right)} + 3 t u{\left(t \right)} \sin{\left(t \right)} - 4 t v{\left(t \right)} \sin{\left(t \right)} \cos{\left(t \right)} - 3 u{\left(t \right)} \cos{\left(t \right)} - 2 v{\left(t \right)} \sin^{2}{\left(t \right)} - v{\left(t \right)} \cos{\left(t \right)}\right)}{24} + \frac{h^{3} \left(- 2 t u{\left(t \right)} \cos{\left(t \right)} - t v{\left(t \right)} \sin^{2}{\left(t \right)} - u{\left(t \right)} \sin{\left(t \right)} - v{\left(t \right)} \sin{\left(t \right)}\right)}{6} + \frac{h^{2} \left(- t u{\left(t \right)} \sin{\left(t \right)} + v{\left(t \right)} \cos{\left(t \right)}\right)}{2} + h v{\left(t \right)} \sin{\left(t \right)} + u{\left(t \right)}$

Could use sympy subs to evaluate these expressions at particular values of t,u,v,h:

In [44]:
u_poly.subs( {u(t):7.3, v(t):-5.1, h:.01, t:2.5} )
Out[44]:
$\displaystyle 7.26914153913721$

But this is very slow.

Question about range(len()):

In [45]:
len([1,2,3,4,5,5,5,5])
Out[45]:
8
In [46]:
mylist = [1,2,3,4,50,50,50,50]
In [47]:
# iteration if we just need the values
for item in mylist:
    print(item)

1
2
3
4
50
50
50
50
In [48]:
# iteration if we need the index as well as the values
for i in range(len(mylist)): 
    print(i,mylist[i])

0 1
1 2
2 3
3 4
4 50
5 50
6 50
7 50

We are done with symbolic computation¶

We have Taylor polynomial approximations to u(t+h),v(t+h) given t,u(t),v(t),h

In [49]:
u_poly
Out[49]:
$\displaystyle \frac{h^{4} \left(t^{2} u{\left(t \right)} \sin^{2}{\left(t \right)} + 3 t u{\left(t \right)} \sin{\left(t \right)} - 4 t v{\left(t \right)} \sin{\left(t \right)} \cos{\left(t \right)} - 3 u{\left(t \right)} \cos{\left(t \right)} - 2 v{\left(t \right)} \sin^{2}{\left(t \right)} - v{\left(t \right)} \cos{\left(t \right)}\right)}{24} + \frac{h^{3} \left(- 2 t u{\left(t \right)} \cos{\left(t \right)} - t v{\left(t \right)} \sin^{2}{\left(t \right)} - u{\left(t \right)} \sin{\left(t \right)} - v{\left(t \right)} \sin{\left(t \right)}\right)}{6} + \frac{h^{2} \left(- t u{\left(t \right)} \sin{\left(t \right)} + v{\left(t \right)} \cos{\left(t \right)}\right)}{2} + h v{\left(t \right)} \sin{\left(t \right)} + u{\left(t \right)}$
In [50]:
v_poly
Out[50]:
$\displaystyle \frac{h^{4} \left(2 t^{2} u{\left(t \right)} \cos{\left(t \right)} + t^{2} v{\left(t \right)} \sin^{2}{\left(t \right)} + 4 t u{\left(t \right)} \sin{\left(t \right)} + t v{\left(t \right)} \sin{\left(t \right)} - 3 v{\left(t \right)} \cos{\left(t \right)}\right)}{24} + \frac{h^{3} \left(t^{2} u{\left(t \right)} \sin{\left(t \right)} - t v{\left(t \right)} \cos{\left(t \right)} - 2 v{\left(t \right)} \sin{\left(t \right)}\right)}{6} + \frac{h^{2} \left(- t v{\left(t \right)} \sin{\left(t \right)} - u{\left(t \right)}\right)}{2} - h t u{\left(t \right)} + v{\left(t \right)}$

For speed we need to move from sympy to numpy

To do this, we can use a sympy function called lambdify which generates a numpy-friendly function from a sympy expression

In [51]:
expr = t**2 + h
expr
Out[51]:
$\displaystyle h + t^{2}$
In [52]:
# arguments of lambdify are the argument(s) of the generated function, 
# the expression(s) to be lambdified, and a keyword "numpy"
foo = sp.lambdify( (t,h), expr, 'numpy' )
In [53]:
foo(5.5, .01)
Out[53]:
30.26
In [54]:
import numpy as np
foo(5.5, np.array([.1,.01,.001]))  # function will work for numpy array inputs
Out[54]:
array([30.35 , 30.26 , 30.251])
In [55]:
# Here we lambdify our Taylor polynomial expressions:

nextuv = sp.lambdify( (t,h,u(t),v(t)), (u_poly,v_poly), 'numpy' )
In [56]:
# Then we can evaluate them at a particular (t,h,u,v) as follows:

nextuv(2.5, .0001, 7.3, -4.2)
Out[56]:
(7.2997486039181, -4.201825005076226)
In [62]:
# Now we just need to use it to generate our approximation to the solution of the differential equation:

t = 0
u = 2
v = 5
tfinal = 3
nsteps = 20
h = (tfinal-t)/nsteps

print(f'{t:5.3f}  {u:13.10f} {v:13.10f}')
for i in range(nsteps):
    u,v = nextuv(t,h,u,v)
    t += h
    print(f'{t:5.3f}  {u:13.10f} {v:13.10f}')

0.000   2.0000000000  5.0000000000
0.150   2.0560179687  4.9771835938
0.300   2.2212440601  4.9049865562
0.450   2.4869270577  4.7724993379
0.600   2.8377303003  4.5626566302
0.750   3.2509683472  4.2540238877
0.900   3.6962148978  3.8234451495
1.050   4.1358288620  3.2495407997
1.200   4.5268908433  2.5168527675
1.350   4.8248172806  1.6201744028
1.500   4.9885238948  0.5683382293
1.650   4.9865081720 -0.6134356322
1.800   4.8027479626 -1.8835007241
1.950   4.4410489248 -3.1866334548
2.100   3.9265760138 -4.4598028624
2.250   3.3038320401 -5.6402091355
2.400   2.6312022288 -6.6743159579
2.550   1.9731212443 -7.5262621320
2.700   1.3916247937 -8.1841381152
2.850   0.9392698602 -8.6631704450
3.000   0.6550699786 -9.0057315860
In [63]:
20//4, 20/4
Out[63]:
(5, 5.0)
In [ ]: