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.



In [2]:

    
import math

math.sqrt(3)









    Out[2]:





1.7320508075688772



In [3]:

    
math.sqrt(8)









    Out[3]:





2.8284271247461903

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



In [4]:

    
import sympy
sympy.sqrt(3)









    Out[4]:





$\displaystyle \sqrt{3}$

SymPy can even simplify symbolic computations



In [5]:

    
sympy.sqrt(8)









    Out[5]:





$\displaystyle 2 \sqrt{2}$



In [6]:

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









    Out[6]:





$\displaystyle x + 2 y$

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



In [7]:

    
expr + 1









    Out[7]:





$\displaystyle x + 2 y + 1$



In [8]:

    
expr - x









    Out[8]:





$\displaystyle 2 y$

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



In [9]:

    
x*expr









    Out[9]:





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



In [10]:

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









    Out[10]:





$\displaystyle x^{2} + 2 x y$



In [11]:

    
factor(expanded_expr)









    Out[11]:





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



In [12]:

    
from sympy import diff, sin, exp

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









    Out[12]:





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



In [13]:

    
from sympy import limit

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









    Out[13]:





$\displaystyle 1$

Exercise

Solve $x^2 - 2 = 0$ using sympy.solve



In [13]:

    
# Type solution here
from sympy import solve

Pretty printing



In [14]:

    
from sympy import init_printing, Integral, sqrt

init_printing(use_latex='mathjax')



In [15]:

    
Integral(sqrt(1/x), x)









    Out[15]:





$\displaystyle \int \sqrt{\frac{1}{x}}\, dx$



In [16]:

    
from sympy import latex

latex(Integral(sqrt(1/x), x))









    Out[16]:





'\\int \\sqrt{\\frac{1}{x}}\\, dx'

More symbols. Exercise: fix the following piece of code



In [18]:

    
# 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                                 Traceback (most recent call last)
<ipython-input-18-88f4dc46acc5> in <module>
----> 1 expr2 = x + 2*y +3*z

NameError: name 'z' is not defined

Exercise

Solve $x + 2*y + 3*z$ for $x$



In [19]:

    
# Solution here
from sympy import solve

Difference between symbol name and python variable name



In [20]:

    
x, y = symbols("y z")



In [21]:

    
x









    Out[21]:





$\displaystyle y$



In [22]:

    
y









    Out[22]:





$\displaystyle z$



In [23]:

    
# NBVAL_SKIP
# The following code will error until the code in cell 16 above is
# fixed
z









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-23-3a710d2a84f8> in <module>
----> 1 z

NameError: name 'z' is not defined

Symbol names can be more than one character long



In [24]:

    
crazy = symbols('unrelated')

crazy + 1









    Out[24]:





$\displaystyle unrelated + 1$



In [25]:

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

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



In [26]:

    
print(expr)

How do we get 3?



In [27]:

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









    Out[27]:





$\displaystyle 3$

Equalities



In [28]:

    
x + 1 == 4









    Out[28]:





False



In [29]:

    
from sympy import Eq

Eq(x + 1, 4)









    Out[29]:





$\displaystyle x + 1 = 4$

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



In [30]:

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









    Out[30]:





False



In [31]:

    
from sympy import simplify

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

simplify(a-b)









    Out[31]:





$\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



In [32]:

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









    Out[32]:





$\displaystyle 36$



In [33]:

    
from sympy import sympify

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









    Out[33]:





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



In [34]:

    
expr.subs(x, 2)









    Out[34]:





$\displaystyle \frac{19}{2}$



In [35]:

    
expr = sqrt(8)



In [36]:

    
expr









    Out[36]:





$\displaystyle 2 \sqrt{2}$



In [37]:

    
expr.evalf()









    Out[37]:





$\displaystyle 2.82842712474619$



In [38]:

    
from sympy import pi

pi.evalf(100)









    Out[38]:





$\displaystyle 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068$



In [39]:

    
from sympy import cos

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









    Out[39]:





$\displaystyle 0.0874989834394464$

Exercise



In [41]:

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









    Out[41]:

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



In [42]:

    
from sympy import Function

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









    Out[42]:





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



In [43]:

    
f(x).diff()









    Out[43]:





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



In [44]:

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









    Out[44]:





$\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)}$



In [45]:

    
from sympy import dsolve

dsolve(diffeq, f(x))









    Out[45]:





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

Finite Differences



In [48]:

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









    Out[48]:





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



In [49]:

    
from sympy import Symbol

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









    Out[49]:





$\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!