Sympy¶

In [1]:

%config InlineBackend.figure_format = 'svg'
import numpy as np
import matplotlib.pyplot as plt

Let us import everything from sympy. We will also enable pretty printing.

In [2]:

from sympy import *
init_printing()

We must define symbols which will be treated symbolically

In [3]:

x = symbols('x')
type(x)

Out[3]:

sympy.core.symbol.Symbol

Let us define an expression

In [4]:

f = x*sin(pi*x) + tan(pi*x)
print(f)

x*sin(pi*x) + tan(pi*x)

Differentiation. We can perform symbolic differentiation.

In [5]:

diff(f,x)

Out[5]:

$\displaystyle \pi x \cos{\left(\pi x \right)} + \pi \left(\tan^{2}{\left(\pi x \right)} + 1\right) + \sin{\left(\pi x \right)}$

Compute second derivative

In [6]:

diff(f,x,2)

Out[6]:

$\displaystyle \pi \left(- \pi x \sin{\left(\pi x \right)} + 2 \pi \left(\tan^{2}{\left(\pi x \right)} + 1\right) \tan{\left(\pi x \right)} + 2 \cos{\left(\pi x \right)}\right)$

Product of two functions. Define a more complicated expression as product of two expressions

$$ h(x) = f(x) g(x) $$

and differentiate it

In [7]:

g = exp(x)
h = f*g
diff(h,x)

Out[7]:

$\displaystyle \left(x \sin{\left(\pi x \right)} + \tan{\left(\pi x \right)}\right) e^{x} + \left(\pi x \cos{\left(\pi x \right)} + \pi \left(\tan^{2}{\left(\pi x \right)} + 1\right) + \sin{\left(\pi x \right)}\right) e^{x}$

Function composition. Composition of two functions

$$ h(x) = (f \circ g)(x) = f(g(x)) $$

Apply this to the case $f(x) = \sin(x)$, $g(x) = \cos(x)$.

In [8]:

g = cos(x)
h = sin(g)
print('h =',h)
diff(h,x)

h = sin(cos(x))

Out[8]:

$\displaystyle - \sin{\left(x \right)} \cos{\left(\cos{\left(x \right)} \right)}$

If the two functions $f,g$ have already been defined and we want to compose, then we can use substitution,

In [9]:

f = sin(x)
g = cos(x)
h = f.subs({'x':g})
print('h =',h)
diff(h,x)

h = sin(cos(x))

Out[9]:

$\displaystyle - \sin{\left(x \right)} \cos{\left(\cos{\left(x \right)} \right)}$

Integration. We can also compute indefinite integrals

In [10]:

print(f)
integrate(f,x)

sin(x)

Out[10]:

$\displaystyle - \cos{\left(x \right)}$

We can get definite integrals.

In [11]:

integrate(f,(x,0,1))

Out[11]:

$\displaystyle 1 - \cos{\left(1 \right)}$

Taylor expansion. Compute Taylor expansion of $f(x)$ around $x=0$

In [12]:

print(f)
series(f,x,0,8)

sin(x)

Out[12]:

$\displaystyle x - \frac{x^{3}}{6} + \frac{x^{5}}{120} - \frac{x^{7}}{5040} + O\left(x^{8}\right)$

Compute Taylor expansion of $h(x)$ around $x=0$

In [13]:

print(h)
series(h,x,0,8)

sin(cos(x))

Out[13]:

$\displaystyle \sin{\left(1 \right)} - \frac{x^{2} \cos{\left(1 \right)}}{2} + x^{4} \left(- \frac{\sin{\left(1 \right)}}{8} + \frac{\cos{\left(1 \right)}}{24}\right) + x^{6} \left(\frac{7 \cos{\left(1 \right)}}{360} + \frac{\sin{\left(1 \right)}}{48}\right) + O\left(x^{8}\right)$

Expression is not a function¶

In [14]:

p = (x + 1)/(x**2 + 2)

We cannot evaluate p at a numerical value. To do that, we can substitute the symbol x with a numerical value.

In [15]:

p.subs({'x':1})

Out[15]:

$\displaystyle \frac{2}{3}$

In [16]:

p.subs({'x':1.0})

Out[16]:

$\displaystyle 0.666666666666667$

The result depends on whether we substitute with an integer or a float. You can substitute with a rational number like this

In [17]:

p.subs({'x':1//2})

Out[17]:

$\displaystyle \frac{1}{2}$

For how to make python functions out of expressions, see below.

Compute roots¶

Let us find roots of quadratic equation

$$ a x^2 + b x + c = 0 $$

In [18]:

x, a, b, c = symbols('x a b c')
eq = a*x**2 + b*x + c
roots(eq,x)

Out[18]:

$\displaystyle \left\{ - \frac{b}{2 a} - \frac{\sqrt{- 4 a c + b^{2}}}{2 a} : 1, \ - \frac{b}{2 a} + \frac{\sqrt{- 4 a c + b^{2}}}{2 a} : 1\right\}$

In [19]:

x, a = symbols('x a')
eq = (x - a)**2
roots(eq,x)

Out[19]:

$\displaystyle \left\{ a : 2\right\}$

Expand, factor, simplify¶

In [20]:

x, a = symbols('x a')
expand((x-a)**2)

Out[20]:

$\displaystyle a^{2} - 2 a x + x^{2}$

In [21]:

factor(x**2 - 2*a*x + a**2)

Out[21]:

$\displaystyle \left(- a + x\right)^{2}$

In [22]:

f = (x + x**2)/(x*sin(x)**2 + x*cos(x)**2)
simplify(f)

Out[22]:

$\displaystyle x + 1$

Get Latex code¶

In [23]:

x = symbols('x')
p = (cos(x)-sin(x))**3
expand(p)

Out[23]:

$\displaystyle - \sin^{3}{\left(x \right)} + 3 \sin^{2}{\left(x \right)} \cos{\left(x \right)} - 3 \sin{\left(x \right)} \cos^{2}{\left(x \right)} + \cos^{3}{\left(x \right)}$

In [24]:

print_latex(expand(p))

- \sin^{3}{\left(x \right)} + 3 \sin^{2}{\left(x \right)} \cos{\left(x \right)} - 3 \sin{\left(x \right)} \cos^{2}{\left(x \right)} + \cos^{3}{\left(x \right)}

Solve a linear system¶

$$ x_1 + x_2 = a, \qquad x_1 - x_2 = b $$

In [25]:

x1, x2, a, b = symbols('x1 x2 a b')
e1 = x1 + x2 - a
e2 = x1 - x2 - b
solve([e1, e2],[x1,x2])

Out[25]:

$\displaystyle \left\{ x_{1} : \frac{a}{2} + \frac{b}{2}, \ x_{2} : \frac{a}{2} - \frac{b}{2}\right\}$

We can capture the result as a dictionary.

In [26]:

r = solve([e1, e2],[x1,x2])
print('x1 = ', r[x1])
print('x2 = ', r[x2])

x1 =  a/2 + b/2
x2 =  a/2 - b/2

Example: create Python function from sympy expression¶

Lets first define a function and get its derivative.

In [27]:

x = symbols('x')
f = x * sin(50*x) * exp(x)
g = diff(f,x)

We cannot evaluate f and g since they are not python functions. We first create Python functions out of the symbolic expressions and then plot them.

In [28]:

ffun = lambdify(x,f)
gfun = lambdify(x,g)

We can now evaluate these functions at some argument

In [29]:

print('f(1) =',ffun(1.0))
print('g(1) =',gfun(1.0))

f(1) = -0.7132087970679899
g(1) = 129.72606342238427

We can now use these to plot graphs of the functions

In [30]:

xx = np.linspace(0.0,1.0,200)
plt.figure(figsize=(9,4))
plt.subplot(121)
plt.plot(xx,ffun(xx))
plt.grid(True)
plt.title('f(x) = '+str(f))
plt.subplot(122)
plt.plot(xx,gfun(xx))
plt.grid(True)
plt.title('Derivative of f(x)');

Example: Truncation error of FD scheme¶

Approximate second derivative using finite difference scheme

$$ f''(x) \approx \frac{f(x-h) - 2 f(x) + f(x+h)}{h^2} $$

We perform Taylor expansion around $x$ to find the error in this approximation: write

$$ f(x \pm h) = f(x) \pm h f'(x) + \frac{h^2}{2} f''(x) \pm \ldots $$

and simplify the rhs. This can be done by sympy.

In [31]:

x, h = symbols("x,h")
f = Function("f")
T = lambda h: (f(x-h) - 2*f(x) + f(x+h))/(h**2)
series(T(h), h, x0=0, n=6)

Out[31]:

$\displaystyle \left. \frac{d^{2}}{d \xi_{1}^{2}} f{\left(\xi_{1} \right)} \right|_{\substack{ \xi_{1}=x }} + \frac{h^{2} \left. \frac{d^{4}}{d \xi_{1}^{4}} f{\left(\xi_{1} \right)} \right|_{\substack{ \xi_{1}=x }}}{12} + \frac{h^{4} \left. \frac{d^{6}}{d \xi_{1}^{6}} f{\left(\xi_{1} \right)} \right|_{\substack{ \xi_{1}=x }}}{360} + O\left(h^{6}\right)$

The above result shows that the FD formula is equal to

$$ \frac{f(x-h) - 2 f(x) + f(x+h)}{h^2} = f''(x) + \frac{h^2}{12} f^{(4)}(x) + \frac{h^4}{360} f^{(6)}(x) + O(h^6) $$

The leading error term is $O(h^2)$

$$ \textrm{error} = \frac{f(x-h) - 2 f(x) + f(x+h)}{h^2} - f''(x) = \frac{h^2}{12} f^{(4)}(x) + \frac{h^4}{360} f^{(6)}(x) + O(h^6) = O(h^2) $$

so that the formula is second order accurate.