import math
math.sqrt(3)1.7320508075688772Symbolic 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.
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 solvePretty 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*zNameError: name 'z' is not definedExercise
Solve \(x + 2*y + 3*z\) for \(x\)
# Solution here
from sympy import solveDifference 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
zNameError: name 'z' is not definedSymbol names can be more than one character long
crazy = symbols('unrelated')
crazy + 1\(\displaystyle unrelated + 1\)
x = symbols("x")
expr = x + 1
x = 2What happens when I print expr now? Does it print 3?
print(expr)x + 1How do we get 3?
x = symbols("x")
expr = x + 1
expr.subs(x, 2)\(\displaystyle 3\)
Equalities
x + 1 == 4Falsefrom 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 + 1Falsefrom 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!