Symbolic Computation

Symbolic computation deals with symbols, representing them exactly, instead of numerical approximations (floating point).

We will start with the following borrowed tutorial to introduce the concepts of SymPy. Devito uses SymPy heavily and builds upon it in its DSL.

import math

math.sqrt(3)
1.7320508075688772
math.sqrt(8)
2.8284271247461903

\(\sqrt(8) = 2\sqrt(2)\), but it’s hard to see that here

import sympy
sympy.sqrt(3)

\(\displaystyle \sqrt{3}\)

SymPy can even simplify symbolic computations

sympy.sqrt(8)

\(\displaystyle 2 \sqrt{2}\)

from sympy import symbols
x, y = symbols('x y')
expr = x + 2*y
expr

\(\displaystyle x + 2 y\)

Note that simply adding two symbols creates an expression. Now let’s play around with it.

expr + 1

\(\displaystyle x + 2 y + 1\)

expr - x

\(\displaystyle 2 y\)

Note that expr - x was not x + 2y -x

x*expr

\(\displaystyle x \left(x + 2 y\right)\)

from sympy import expand, factor
expanded_expr = expand(x*expr)
expanded_expr

\(\displaystyle x^{2} + 2 x y\)

factor(expanded_expr)

\(\displaystyle x \left(x + 2 y\right)\)

from sympy import diff, sin, exp

diff(sin(x)*exp(x), x)

\(\displaystyle e^{x} \sin{\left(x \right)} + e^{x} \cos{\left(x \right)}\)

from sympy import limit

limit(sin(x)/x, x, 0)

\(\displaystyle 1\)

Exercise

Solve \(x^2 - 2 = 0\) using sympy.solve

# Type solution here
from sympy import solve

Pretty printing

from sympy import init_printing, Integral, sqrt

init_printing(use_latex='mathjax')
Integral(sqrt(1/x), x)

\(\displaystyle \int \sqrt{\frac{1}{x}}\, dx\)

from sympy import latex

latex(Integral(sqrt(1/x), x))
'\\int \\sqrt{\\frac{1}{x}}\\, dx'

More symbols. Exercise: fix the following piece of code

# NBVAL_SKIP
# The following piece of code is supposed to fail as it is
# The exercise is to fix the code

expr2 = x + 2*y +3*z
NameError: name 'z' is not defined

Exercise

Solve \(x + 2*y + 3*z\) for \(x\)

# Solution here
from sympy import solve

Difference between symbol name and python variable name

x, y = symbols("y z")
x

\(\displaystyle y\)

y

\(\displaystyle z\)

# NBVAL_SKIP
# The following code will error until the code in cell 16 above is
# fixed
z
NameError: name 'z' is not defined

Symbol names can be more than one character long

crazy = symbols('unrelated')

crazy + 1

\(\displaystyle unrelated + 1\)

x = symbols("x")
expr = x + 1
x = 2

What happens when I print expr now? Does it print 3?

print(expr)
x + 1

How do we get 3?

x = symbols("x")
expr = x + 1
expr.subs(x, 2)

\(\displaystyle 3\)

Equalities

x + 1 == 4
False
from sympy import Eq

Eq(x + 1, 4)

\(\displaystyle x + 1 = 4\)

Suppose we want to ask whether \((x + 1)^2 = x^2 + 2x + 1\)

(x + 1)**2 == x**2 + 2*x + 1
False
from sympy import simplify

a = (x + 1)**2
b = x**2 + 2*x + 1

simplify(a-b)

\(\displaystyle 0\)

Exercise

Write a function that takes two expressions as input, and returns a tuple of two booleans. The first if they are equal symbolically, and the second if they are equal mathematically.

More operations

z = symbols("z")
expr = x**3 + 4*x*y - z
expr.subs([(x, 2), (y, 4), (z, 0)])

\(\displaystyle 36\)

from sympy import sympify

str_expr = "x**2 + 3*x - 1/2"
expr = sympify(str_expr)
expr

\(\displaystyle x^{2} + 3 x - \frac{1}{2}\)

expr.subs(x, 2)

\(\displaystyle \frac{19}{2}\)

expr = sqrt(8)
expr

\(\displaystyle 2 \sqrt{2}\)

expr.evalf()

\(\displaystyle 2.82842712474619\)

from sympy import pi

pi.evalf(100)

\(\displaystyle 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068\)

from sympy import cos

expr = cos(2*x)
expr.evalf(subs={x: 2.4})

\(\displaystyle 0.0874989834394464\)

Exercise

from IPython.core.display import Image 
Image(filename='figures/comic.png')

Write a function that takes a symbolic expression (like pi), and determines the first place where 789 appears. Tip: Use the string representation of the number. Python starts counting at 0, but the decimal point offsets this

Solving an ODE

from sympy import Function

f, g = symbols('f g', cls=Function)
f(x)

\(\displaystyle f{\left(x \right)}\)

f(x).diff()

\(\displaystyle \frac{d}{d x} f{\left(x \right)}\)

diffeq = Eq(f(x).diff(x, x) - 2*f(x).diff(x) + f(x), sin(x))
diffeq

\(\displaystyle f{\left(x \right)} - 2 \frac{d}{d x} f{\left(x \right)} + \frac{d^{2}}{d x^{2}} f{\left(x \right)} = \sin{\left(x \right)}\)

from sympy import dsolve

dsolve(diffeq, f(x))

\(\displaystyle f{\left(x \right)} = \left(C_{1} + C_{2} x\right) e^{x} + \frac{\cos{\left(x \right)}}{2}\)

Finite Differences


f = Function('f')
dfdx = f(x).diff(x)
dfdx.as_finite_difference()

\(\displaystyle - f{\left(x - \frac{1}{2} \right)} + f{\left(x + \frac{1}{2} \right)}\)

from sympy import Symbol

d2fdx2 = f(x).diff(x, 2)
h = Symbol('h')
d2fdx2.as_finite_difference(h)

\(\displaystyle - \frac{2 f{\left(x \right)}}{h^{2}} + \frac{f{\left(- h + x \right)}}{h^{2}} + \frac{f{\left(h + x \right)}}{h^{2}}\)

Now that we have seen some relevant features of vanilla SymPy, let’s move on to Devito, which could be seen as SymPy finite differences on steroids!

Back to top