SymPy

Note

You can download this example as a Python script: sympy.py or Jupyter notebook: sympy.ipynb.

Learning Objectives

After completing this chapter readers will be able to:

• Write mathematical expressions with symbols and functions using SymPy.

• Print different forms of expressions and equations with SymPy.

• Differentiate mathematical expressions using SymPy.

• Evaluate mathematical expressions using SymPy.

• Create matrices and do linear algebra using SymPy.

• Solve a linear system of equations with SymPy.

• Simplify mathematical expressions with SymPy.

Introduction

SymPy is an open source, collaboratively developed computer algebra system (CAS) written in Python. It will be used extensively for manipulating symbolic expressions and equations. All of the mathematics needed to formulate the equations of motion of multibody systems can be done with pencil and paper, but the bookkeeping becomes extremely tedious and error prone for systems with even a small number of bodies. SymPy lets a computer handle the tedious aspects (e.g. differentiation or solving linear systems of equations) and reduces the errors one would encounter with pencil and paper. This chapter introduces SymPy and the primary SymPy features we will be using.

Import and Setup

I will import SymPy as follows throughout this book:

import sympy as sm

Since SymPy works with mathematical symbols it’s nice to view SymPy objects in a format that is similar to the math in a textbook. Executing init_printing() at the beginning of your Jupyter Notebook will ensure that SymPy objects render as typeset mathematics. I use the use_latex='mathjax' argument here to disable math png image generation, but that keyword argument is not necessary.

sm.init_printing(use_latex='mathjax')

Symbols

Mathematical symbols are created with the symbols() function. A symbol \(a\) is created like so:

a = sm.symbols('a')
a

\[\displaystyle a\]

This symbol object is of the Symbol type:

type(a)

sympy.core.symbol.Symbol

Multiple symbols can be created with one call to symbols() and SymPy recognizes common Greek symbols by their spelled-out name.

b, t, omega, Omega = sm.symbols('b, t, omega, Omega')
b, t, omega, Omega

\[\displaystyle \left( b, \ t, \ \omega, \ \Omega\right)\]

Note that the argument provided to symbols() does not need to match the Python variable name it is assigned to. Using more verbose Python variable names may make code easier to read and understand, especially if there are many mathematical variables that you need to keep track of. Note that the subscripts are recognized too.

pivot_angle, w2 = sm.symbols('alpha1, omega2')
pivot_angle, w2

\[\displaystyle \left( \alpha_{1}, \ \omega_{2}\right)\]

Exercise

Review the SymPy documentation and create symbols \(q_1, q_2, \ldots, q_{10}\) with a very succint call to symbols().

Solution

sm.symbols('q1:11')

\[\displaystyle \left( q_{1}, \ q_{2}, \ q_{3}, \ q_{4}, \ q_{5}, \ q_{6}, \ q_{7}, \ q_{8}, \ q_{9}, \ q_{10}\right)\]

Undefined Functions

You will also work with undefined mathematical functions in addition to symbols. These will play an important role in setting up differential equations, where you typically don’t know the function, but only its derivative(s). You can create arbitrary functions of variables. In this case, you make a function of \(t\). First create the function name:

f = sm.Function('f')
f

f

This is of a type sympy.core.function.UndefinedFunction.

type(f)

sympy.core.function.UndefinedFunction

Now you can create functions of one or more variables like so:

f(t)

\[\displaystyle f{\left(t \right)}\]

Warning

Due to SymPy’s internal implementations, the type of a function with its argument is not defined as expected:

type(f(t))

f

This can be confusing if you are checking types.

The same UndefinedFunction can be used to create multivariate functions:

f(a, b, omega, t)

\[\displaystyle f{\left(a,b,\omega,t \right)}\]

Exercise

Create a function \(H(x, y, z)\).

Solution

x, y, z = sm.symbols('x, y, z')
sm.Function('H')(x, y, z)

\[\displaystyle H{\left(x,y,z \right)}\]

Symbolic Expressions

Now that you have mathematical variables and functions available, they can be used to construct mathematical expressions. The most basic way to construct expressions is with the standard Python operators +, -, *, /, and **. For example:

expr1 = a + b/omega**2
expr1

\[\displaystyle a + \frac{b}{\omega^{2}}\]

An expression will have the type Add, Mul, or Pow:

type(expr1)

sympy.core.add.Add

This is because SymPy stores expressions behind the scenes as a tree. You can inspect this internal representation by using the srepr() function:

sm.srepr(expr1)

"Add(Symbol('a'), Mul(Symbol('b'), Pow(Symbol('omega'), Integer(-2))))"

This is a visual representation of the tree:

This representation is SymPy’s “true” representation of the symbolic expression. SymPy can display this expression in many other representations, for example the typeset mathematical expression you have already seen is one of those representations. This is important to know, because sometimes the expressions are displayed to you in a way that may be confusing and checking the srepr() version can help clear up misunderstandings. See the manipulation section of the SymPy tutorial for more information on this.

Undefined functions can also be used in expressions just like symbols:

expr2 = f(t) + a*omega
expr2

\[\displaystyle a \omega + f{\left(t \right)}\]

SymPy has a large number of elementary and special functions. See the SymPy documentation on functions for more information. For example, here is an expression that uses sin, Abs, and sqrt():

expr3 = a*sm.sin(omega) + sm.Abs(f(t))/sm.sqrt(b)
expr3

\[\displaystyle a \sin{\left(\omega \right)} + \frac{\left|{f{\left(t \right)}}\right|}{\sqrt{b}}\]

Note that Python integers and floats can also be used when constructing expressions:

expr4 = 5*sm.sin(12) + sm.Abs(-1001)/sm.sqrt(89.2)
expr4

\[\displaystyle 5 \sin{\left(12 \right)} + 105.986768359379\]

Warning

Be careful with numbers, as SymPy may not intepret them as expected. For example:

1/2*a

\[\displaystyle 0.5 a\]

Python does the division before it is multiplied by a, thus a floating point value is created. To fix this you can use the S() function to “sympify” numbers:

sm.S(1)/2*a

\[\displaystyle \frac{a}{2}\]

Or you can ensure the symbol comes first in the division operation:

a/2

\[\displaystyle \frac{a}{2}\]

Lastly, an expression of t:

expr5 = t*sm.sin(omega*f(t)) + f(t)/sm.sqrt(t)
expr5

\[\displaystyle t \sin{\left(\omega f{\left(t \right)} \right)} + \frac{f{\left(t \right)}}{\sqrt{t}}\]

Exercise

Create an expression for the normal distribution function:

(1)\[\frac{1}{\sqrt{2\pi\sigma}}e^{\frac{(x-\mu)^2}{2\sigma^2}}\]

Solution

x, s, m = sm.symbols('x, sigma, mu')
sm.exp((x-m)**2/2/s**2)/sm.sqrt(2*sm.pi*s)

\[\displaystyle \frac{\sqrt{2} e^{\frac{\left(- \mu + x\right)^{2}}{2 \sigma^{2}}}}{2 \sqrt{\pi} \sqrt{\sigma}}\]

Notice that SymPy does some minor manipulation of the expression, but it is equivalent to the form shown in the prompt.

Printing

I introduced the srepr() form of SymPy expressions above and mentioned that expressions can have different representations. For the following srepr() form:

sm.srepr(expr3)

"Add(Mul(Symbol('a'), sin(Symbol('omega'))), Mul(Pow(Symbol('b'), Rational(-1, 2)), Abs(Function('f')(Symbol('t')))))"

There is also a standard representation accessed with the repr() function:

repr(expr3)

'a*sin(omega) + Abs(f(t))/sqrt(b)'

This form matches what you typically would type to create the expression and it returns a string. The print() function will display that string:

print(expr3)

a*sin(omega) + Abs(f(t))/sqrt(b)

SymPy also has a “pretty printer” (pprint()) that makes use of unicode symbols to provide a form that more closely resembles typeset math:

sm.pprint(expr3)

           │f(t)│
a⋅sin(ω) + ──────
             √b  

Lastly, the following lines show how SymPy expressions can be represented as LaTeX code using sympy.printing.latex.latex(). The double backslashes are present because double backslashes represent the escape character in Python strings.

sm.latex(expr3)

'a \\sin{\\left(\\omega \\right)} + \\frac{\\left|{f{\\left(t \\right)}}\\right|}{\\sqrt{b}}'

print(sm.latex(expr3))

a \sin{\left(\omega \right)} + \frac{\left|{f{\left(t \right)}}\right|}{\sqrt{b}}

Warning

When you are working with long expressions, which will be the case in this course, there is no need to print them to the screen. In fact, printing them to the screen make take a long time and fill your entire notebook with an unreadable mess.

Exercise

Print the normal distribution expression

(2)\[\frac{1}{\sqrt{2\pi\sigma}}e^{\frac{(x-\mu)^2}{2\sigma^2}}\]

as a LaTeX string inside an equation environment.

Solution

x, s, m = sm.symbols('x, sigma, mu')
print(sm.latex(sm.exp((x-m)**2/2/s**2)/sm.sqrt(2*sm.pi*s),
               mode='equation'))

\begin{equation}\frac{\sqrt{2} e^{\frac{\left(- \mu + x\right)^{2}}{2 \sigma^{2}}}}{2 \sqrt{\pi} \sqrt{\sigma}}\end{equation}

Differentiating

One of the most tedious tasks in formulating equations of motion is differentiating complex trigonometric expressions. SymPy can calculate derivatives effortlessly. The diff() SymPy function takes an undefined function or an expression and differentiates it with respect to the symbol provided as the second argument:

sm.diff(f(t), t)

\[\displaystyle \frac{d}{d t} f{\left(t \right)}\]

All functions and expressions also have a .diff() method which can be used like so (many SymPy functions exist as standalone functions and methods):

f(t).diff(t)

\[\displaystyle \frac{d}{d t} f{\left(t \right)}\]

expr3 is a more complicated expression:

expr3

\[\displaystyle a \sin{\left(\omega \right)} + \frac{\left|{f{\left(t \right)}}\right|}{\sqrt{b}}\]

It can be differentiated, for example, with respect to \(b\):

expr3.diff(b)

\[\displaystyle - \frac{\left|{f{\left(t \right)}}\right|}{2 b^{\frac{3}{2}}}\]

You can also calculate partial derivatives with respect to successive variables. If you want to first differentiate with respect to \(b\) and then with respect to \(t\) as in the following operation:

(3)\[\frac{\partial^2 h(a, \omega, t, b)}{\partial t \partial b}\]

where:

(4)\[h(a, \omega, t, b) = \displaystyle a \sin{\left(\omega \right)} + \frac{\left|{f{\left(t \right)}}\right|}{\sqrt{b}}\]

then you can use successive arguments to .diff():

expr3.diff(b, t)

\[\displaystyle - \frac{\left(\operatorname{re}{\left(f{\left(t \right)}\right)} \frac{d}{d t} \operatorname{re}{\left(f{\left(t \right)}\right)} + \operatorname{im}{\left(f{\left(t \right)}\right)} \frac{d}{d t} \operatorname{im}{\left(f{\left(t \right)}\right)}\right) \operatorname{sign}{\left(f{\left(t \right)} \right)}}{2 b^{\frac{3}{2}} f{\left(t \right)}}\]

Note that the answer includes real and imaginary components and the signum function.

Warning

SymPy assumes all symbols are complex-valued unless told otherwise. You can attach assumptions to symbols to force them to be real, positive, negative, etc. For example, compare these three outputs:

h = sm.Function('h')
sm.Abs(h(t)).diff(t)

\[\displaystyle \frac{\left(\operatorname{re}{\left(h{\left(t \right)}\right)} \frac{d}{d t} \operatorname{re}{\left(h{\left(t \right)}\right)} + \operatorname{im}{\left(h{\left(t \right)}\right)} \frac{d}{d t} \operatorname{im}{\left(h{\left(t \right)}\right)}\right) \operatorname{sign}{\left(h{\left(t \right)} \right)}}{h{\left(t \right)}}\]

h = sm.Function('h', real=True)
sm.Abs(h(t)).diff(t)

\[\displaystyle \operatorname{sign}{\left(h{\left(t \right)} \right)} \frac{d}{d t} h{\left(t \right)}\]

h = sm.Function('h', real=True, positive=True)
sm.Abs(h(t)).diff(t)

\[\displaystyle \frac{d}{d t} h{\left(t \right)}\]

Sometimes you may need to add assumptions to variables, but in general it will not be necessary. Read more about assumptions in SymPy’s guide.

Exercise

Differentiate expr5 above using this operator:

(5)\[\frac{\partial^2}{\partial \omega \partial t}\]

Solution

First show expr5:

expr5

\[\displaystyle t \sin{\left(\omega f{\left(t \right)} \right)} + \frac{f{\left(t \right)}}{\sqrt{t}}\]

The twice partial derivative is:

expr5.diff(t, omega)

\[\displaystyle - \omega t f{\left(t \right)} \sin{\left(\omega f{\left(t \right)} \right)} \frac{d}{d t} f{\left(t \right)} + t \cos{\left(\omega f{\left(t \right)} \right)} \frac{d}{d t} f{\left(t \right)} + f{\left(t \right)} \cos{\left(\omega f{\left(t \right)} \right)}\]

or you can chain .diff() calls:

expr5.diff(t).diff(omega)

\[\displaystyle - \omega t f{\left(t \right)} \sin{\left(\omega f{\left(t \right)} \right)} \frac{d}{d t} f{\left(t \right)} + t \cos{\left(\omega f{\left(t \right)} \right)} \frac{d}{d t} f{\left(t \right)} + f{\left(t \right)} \cos{\left(\omega f{\left(t \right)} \right)}\]

Evaluating Symbolic Expressions

SymPy expressions can be evaluated numerically in several ways. The xreplace() method allows substitution of exact symbols or sub-expressions. First create a dictionary that maps symbols, functions or sub-expressions to the replacements:

repl = {omega: sm.pi/4, a: 2, f(t): -12, b: 25}

This dictionary can then be passed to .xreplace():

expr3.xreplace(repl)

\[\displaystyle \sqrt{2} + \frac{12}{5}\]

Notice how the square root and fraction do not automatically reduce to their decimal equivalents. To do so, you must use the evalf() method. This method will evaluate an expression to an arbitrary number of decimal points. You provide the number of decimal places and the substitution dictionary to evaluate:

expr3.evalf(n=10, subs=repl)

\[\displaystyle 3.814213562\]

type(expr3.evalf(n=10, subs=repl))

sympy.core.numbers.Float

Note that this is a SymPy Float object, which is a special object that can have an arbitrary number of decimal places, for example here is the expression evaluated to 80 decimal places:

expr3.evalf(n=80, subs=repl)

\[\displaystyle 3.814213562373095048801688724209698078569671875376948073176679737990732478462107\]

To convert this to Python floating point number, use float():

float(expr3.evalf(n=300, subs=repl))

\[\displaystyle 3.81421356237309\]

type(float(expr3.evalf(n=300, subs=repl)))

float

This value is a machine precision floating point value and can be used with standard Python functions that operate on floating point numbers.

To obtain machine precision floating point numbers directly and with more flexibility, it is better to use the lambdify() function to convert the expression to a Python function. When using lambdify(), all symbols and functions should be converted to numbers so first identify what symbols and functions make up the expression.

expr3

\[\displaystyle a \sin{\left(\omega \right)} + \frac{\left|{f{\left(t \right)}}\right|}{\sqrt{b}}\]

\(\omega, a, f(t)\), and \(b\) are all present in the expression. The first argument of lambdify() should be a sequence of all these symbols and functions and the second argument should be the expression.

eval_expr3 = sm.lambdify((omega, a, f(t), b), expr3)

lambdify() generates a Python function and, in this case, we store that function in the variable eval_expr3. You can see what the inputs and outputs of the function are with help():

help(eval_expr3)

Help on function _lambdifygenerated:

_lambdifygenerated(omega, a, _Dummy_38, b)
    Created with lambdify. Signature:
    
    func(omega, a, f, b)
    
    Expression:
    
    a*sin(omega) + Abs(f(t))/sqrt(b)
    
    Source code:
    
    def _lambdifygenerated(omega, a, _Dummy_38, b):
        return a*sin(omega) + abs(_Dummy_38)/sqrt(b)
    
    
    Imported modules:

This function operates on and returns floating point values, for example:

eval_expr3(3.14/4, 2, -12, 25)

\[\displaystyle 3.81365036221073\]

The type of lambdify’s return values will be NumPy floats.

type(eval_expr3(3.14/4, 2, -12, 25))

numpy.float64

These floats are interoperable with Python floats for single values (unlike SymPy Floats) but also support arrays of floats. For example:

eval_expr3(3.14/4, 2, -12, [25, 26, 27])

array([3.81365036, 3.76704398, 3.72305144])

type(eval_expr3(3.14/4, 2, -12, [25, 26, 27]))

numpy.ndarray

More on NumPy arrays of floats will be introduced in a later chapter.

Warning

Python and NumPy floats can be mixed, but avoid mixing SymPy Floats with either.

Note

This distinction between SymPy Float objects and regular Python and NumPy float objects is important. In this case, the Python float and the NumPy float are equivalent. The later will compute much faster because arbitrary precision is not required. In this book, you will almost always want to convert SymPy expressions to machine precision floating point numbers, so use lambdify().

Exercise

Create a symbolic expression representing Newton’s Law of Universal Gravitation. Use lambdify() to evaluate the expression for two mass of 5.972E24 kg and 80 kg at a distance of 6371 km apart to find the gravitational force in Newtons.

Solution

G, m1, m2, r = sm.symbols('G, m1, m2, r')
F = G*m1*m2/r**2
eval_F = sm.lambdify((G, m1, m2, r), F)
eval_F(6.67430E-11, 5.972E24, 80, 6371E3)

\[\displaystyle 785.597874097975\]

Matrices

SymPy supports matrices of expressions and linear algebra. Many of the operations needed in multibody dynamics are more succinctly formulated with matrices and linear algebra. Matrices can be created by passing nested lists to the Matrix() object. For example:

mat1 = sm.Matrix([[a, 2*a], [b/omega, f(t)]])
mat1

\[\begin{split}\displaystyle \left[\begin{matrix}a & 2 a\\\frac{b}{\omega} & f{\left(t \right)}\end{matrix}\right]\end{split}\]

mat2 = sm.Matrix([[1, 2], [3, 4]])
mat2

\[\begin{split}\displaystyle \left[\begin{matrix}1 & 2\\3 & 4\end{matrix}\right]\end{split}\]

All matrices are two dimensional and the number of rows and columns, in that order, are stored in the .shape attribute.

mat1.shape

\[\displaystyle \left( 2, \ 2\right)\]

Individual elements of the matrix can be extracted with the bracket notation taking the row and column indices (remember Python indexes from 0):

mat1[0, 1]

\[\displaystyle 2 a\]

The slice notation can extract rows or columns:

mat1[0, 0:2]

\[\displaystyle \left[\begin{matrix}a & 2 a\end{matrix}\right]\]

mat1[0:2, 1]

\[\begin{split}\displaystyle \left[\begin{matrix}2 a\\f{\left(t \right)}\end{matrix}\right]\end{split}\]

Matrix algebra can be performed. Matrices can be added:

mat1 + mat2

\[\begin{split}\displaystyle \left[\begin{matrix}a + 1 & 2 a + 2\\\frac{b}{\omega} + 3 & f{\left(t \right)} + 4\end{matrix}\right]\end{split}\]

Both the * and the @ operator perform matrix multiplication:

mat1*mat2

\[\begin{split}\displaystyle \left[\begin{matrix}7 a & 10 a\\\frac{b}{\omega} + 3 f{\left(t \right)} & \frac{2 b}{\omega} + 4 f{\left(t \right)}\end{matrix}\right]\end{split}\]

mat1@mat2

\[\begin{split}\displaystyle \left[\begin{matrix}7 a & 10 a\\\frac{b}{\omega} + 3 f{\left(t \right)} & \frac{2 b}{\omega} + 4 f{\left(t \right)}\end{matrix}\right]\end{split}\]

Element-by-element multiplication requires the sympy.hadamard_product() function:

sm.hadamard_product(mat1, mat2)

\[\begin{split}\displaystyle \left[\begin{matrix}a & 4 a\\\frac{3 b}{\omega} & 4 f{\left(t \right)}\end{matrix}\right]\end{split}\]

Note that NumPy uses * for element-by-element multiplication and @ for matrix multiplication, so to avoid possible confusion, use @ for SymPy matrix multiplication.

Differentiation operates on each element of the matrix:

mat3 = sm.Matrix([expr1, expr2, expr3, expr4, expr5])
mat3

\[\begin{split}\displaystyle \left[\begin{matrix}a + \frac{b}{\omega^{2}}\\a \omega + f{\left(t \right)}\\a \sin{\left(\omega \right)} + \frac{\left|{f{\left(t \right)}}\right|}{\sqrt{b}}\\5 \sin{\left(12 \right)} + 105.986768359379\\t \sin{\left(\omega f{\left(t \right)} \right)} + \frac{f{\left(t \right)}}{\sqrt{t}}\end{matrix}\right]\end{split}\]

mat3.diff(a)

\[\begin{split}\displaystyle \left[\begin{matrix}1\\\omega\\\sin{\left(\omega \right)}\\0\\0\end{matrix}\right]\end{split}\]

mat3.diff(t)

\[\begin{split}\displaystyle \left[\begin{matrix}0\\\frac{d}{d t} f{\left(t \right)}\\\frac{\left(\operatorname{re}{\left(f{\left(t \right)}\right)} \frac{d}{d t} \operatorname{re}{\left(f{\left(t \right)}\right)} + \operatorname{im}{\left(f{\left(t \right)}\right)} \frac{d}{d t} \operatorname{im}{\left(f{\left(t \right)}\right)}\right) \operatorname{sign}{\left(f{\left(t \right)} \right)}}{\sqrt{b} f{\left(t \right)}}\\0\\\omega t \cos{\left(\omega f{\left(t \right)} \right)} \frac{d}{d t} f{\left(t \right)} + \sin{\left(\omega f{\left(t \right)} \right)} + \frac{\frac{d}{d t} f{\left(t \right)}}{\sqrt{t}} - \frac{f{\left(t \right)}}{2 t^{\frac{3}{2}}}\end{matrix}\right]\end{split}\]

If you have column vectors \(\bar{v}\) and \(\bar{u}\), the \((i,j)\) entries of the Jacobian of \(\bar{v}\) with respect to the entries in vector \(\bar{u}\) are found with \(\mathbf{J}_{ij} = \frac{\partial v_i}{\partial u_j}\). The Jacobian matrix of vector (column matrix) can be formed with the jacobian() method. This calculates the partial derivatives of each element in the vector with respect to a vector (or sequence) of variables.

mat4 = sm.Matrix([a, b, omega, t])
mat4

\[\begin{split}\displaystyle \left[\begin{matrix}a\\b\\\omega\\t\end{matrix}\right]\end{split}\]

mat3.jacobian(mat4)

\[\begin{split}\displaystyle \left[\begin{matrix}1 & \frac{1}{\omega^{2}} & - \frac{2 b}{\omega^{3}} & 0\\\omega & 0 & a & \frac{d}{d t} f{\left(t \right)}\\\sin{\left(\omega \right)} & - \frac{\left|{f{\left(t \right)}}\right|}{2 b^{\frac{3}{2}}} & a \cos{\left(\omega \right)} & \frac{\left(\operatorname{re}{\left(f{\left(t \right)}\right)} \frac{d}{d t} \operatorname{re}{\left(f{\left(t \right)}\right)} + \operatorname{im}{\left(f{\left(t \right)}\right)} \frac{d}{d t} \operatorname{im}{\left(f{\left(t \right)}\right)}\right) \operatorname{sign}{\left(f{\left(t \right)} \right)}}{\sqrt{b} f{\left(t \right)}}\\0 & 0 & 0 & 0\\0 & 0 & t f{\left(t \right)} \cos{\left(\omega f{\left(t \right)} \right)} & \omega t \cos{\left(\omega f{\left(t \right)} \right)} \frac{d}{d t} f{\left(t \right)} + \sin{\left(\omega f{\left(t \right)} \right)} + \frac{\frac{d}{d t} f{\left(t \right)}}{\sqrt{t}} - \frac{f{\left(t \right)}}{2 t^{\frac{3}{2}}}\end{matrix}\right]\end{split}\]

Exercise

Write your own function that produces a Jacobian given a column matrix of expressions. It should look like:

def jacobian(v, x):
    """Returns the Jacobian of the vector function v with respect to the
    vector of variables x."""
    # fill in your code here
    return J_v_x

Show that it gives the same solution as the above .jacobian() method. Do not use the .jacobian() method in your function.

Solution

def jacobian(v, x):
    """Returns the Jacobian of the vector function v with respect to the
    vector of variables x."""
    diffs = []
    for expr in v:
      for var in x:
         diffs.append(expr.diff(var))
    J_v_x = sm.Matrix(diffs).reshape(len(v), len(x))
    return J_v_x

jacobian(mat3, mat4)

\[\begin{split}\displaystyle \left[\begin{matrix}1 & \frac{1}{\omega^{2}} & - \frac{2 b}{\omega^{3}} & 0\\\omega & 0 & a & \frac{d}{d t} f{\left(t \right)}\\\sin{\left(\omega \right)} & - \frac{\left|{f{\left(t \right)}}\right|}{2 b^{\frac{3}{2}}} & a \cos{\left(\omega \right)} & \frac{\left(\operatorname{re}{\left(f{\left(t \right)}\right)} \frac{d}{d t} \operatorname{re}{\left(f{\left(t \right)}\right)} + \operatorname{im}{\left(f{\left(t \right)}\right)} \frac{d}{d t} \operatorname{im}{\left(f{\left(t \right)}\right)}\right) \operatorname{sign}{\left(f{\left(t \right)} \right)}}{\sqrt{b} f{\left(t \right)}}\\0 & 0 & 0 & 0\\0 & 0 & t f{\left(t \right)} \cos{\left(\omega f{\left(t \right)} \right)} & \omega t \cos{\left(\omega f{\left(t \right)} \right)} \frac{d}{d t} f{\left(t \right)} + \sin{\left(\omega f{\left(t \right)} \right)} + \frac{\frac{d}{d t} f{\left(t \right)}}{\sqrt{t}} - \frac{f{\left(t \right)}}{2 t^{\frac{3}{2}}}\end{matrix}\right]\end{split}\]

Solving Linear Systems

You’ll need to solve linear systems of equations often in this book. SymPy offers a number of ways to do this, but the best way to do so if you know a set of equations are linear in specific variables is the method described below. First, you should confirm you have equations of this form:

(6)\[\begin{split}a_{11} x_1 + a_{12} x_2 + \ldots + a_{1n} x_n + b_1 = 0 \\ a_{21} x_1 + a_{22} x_2 + \ldots + a_{2n} x_n + b_2 = 0 \\ \vdots \\ a_{n1} x_1 + a_{n2} x_2 + \ldots + a_{nn} x_n + b_n = 0\end{split}\]

These equations can be put into matrix form:

(7)\[\mathbf{A}\bar{x} = \bar{b}\]

where:

(8)\[\begin{split}\mathbf{A} = \begin{bmatrix} a_{11} & a_{12} & \ldots & a_{1n} \\ a_{21} & a_{22} & \ldots & a_{2n} \\ \ldots & \ldots & \ldots & \ldots \\ a_{n1} & a_{n2} & \ldots & a_{nn} \end{bmatrix}, \bar{x} = \begin{bmatrix} x_1 \\ x_2 \\ \ldots \\ x_n \end{bmatrix}, \bar{b} = \begin{bmatrix} -b_1 \\ -b_2 \\ \ldots \\ -b_n \end{bmatrix}\end{split}\]

\(\bar{x}\), the solution, is found with matrix inversion (if the matrix is invertible):

(9)\[\bar{x} = \mathbf{A}^{-1}\bar{b}\]

Taking the inverse is not computationally efficient and potentially numerically inaccurate, so some form of Gaussian elmination should be used to solve the system.

To solve with SymPy, start with a column matrix of linear expressions:

a1, a2 = sm.symbols('a1, a2')

exprs = sm.Matrix([
    [a1*sm.sin(f(t))*sm.cos(2*f(t)) + a2 + omega/sm.log(f(t), t) + 100],
    [a1*omega**2 + f(t)*a2 + omega + f(t)**3],
])
exprs

\[\begin{split}\displaystyle \left[\begin{matrix}a_{1} \sin{\left(f{\left(t \right)} \right)} \cos{\left(2 f{\left(t \right)} \right)} + a_{2} + \frac{\omega \log{\left(t \right)}}{\log{\left(f{\left(t \right)} \right)}} + 100\\a_{1} \omega^{2} + a_{2} f{\left(t \right)} + \omega + f^{3}{\left(t \right)}\end{matrix}\right]\end{split}\]

Since we know these two expressions are linear in the \(a_1\) and \(a_2\) variables, the partial derivatives with respect to those two variables will return the linear coefficients. The \(\mathbf{A}\) matrix can be formed in one step with the .jacobian() method:

A = exprs.jacobian([a1, a2])
A

\[\begin{split}\displaystyle \left[\begin{matrix}\sin{\left(f{\left(t \right)} \right)} \cos{\left(2 f{\left(t \right)} \right)} & 1\\\omega^{2} & f{\left(t \right)}\end{matrix}\right]\end{split}\]

The \(\bar{b}\) vector can be formed by setting \(a_1=a_2=0\), leaving the terms that are not linear in \(a_1\) and \(a_2\).

b = -exprs.xreplace({a1: 0, a2: 0})
b

\[\begin{split}\displaystyle \left[\begin{matrix}- \frac{\omega \log{\left(t \right)}}{\log{\left(f{\left(t \right)} \right)}} - 100\\- \omega - f^{3}{\left(t \right)}\end{matrix}\right]\end{split}\]

The inv() method can compute the inverse of A to find the solution:

A.inv() @ b

\[\begin{split}\displaystyle \left[\begin{matrix}- \frac{- \omega - f^{3}{\left(t \right)}}{- \omega^{2} + f{\left(t \right)} \sin{\left(f{\left(t \right)} \right)} \cos{\left(2 f{\left(t \right)} \right)}} + \frac{\left(- \frac{\omega \log{\left(t \right)}}{\log{\left(f{\left(t \right)} \right)}} - 100\right) f{\left(t \right)}}{- \omega^{2} + f{\left(t \right)} \sin{\left(f{\left(t \right)} \right)} \cos{\left(2 f{\left(t \right)} \right)}}\\- \frac{\omega^{2} \left(- \frac{\omega \log{\left(t \right)}}{\log{\left(f{\left(t \right)} \right)}} - 100\right)}{- \omega^{2} + f{\left(t \right)} \sin{\left(f{\left(t \right)} \right)} \cos{\left(2 f{\left(t \right)} \right)}} + \frac{\left(- \omega - f^{3}{\left(t \right)}\right) \sin{\left(f{\left(t \right)} \right)} \cos{\left(2 f{\left(t \right)} \right)}}{- \omega^{2} + f{\left(t \right)} \sin{\left(f{\left(t \right)} \right)} \cos{\left(2 f{\left(t \right)} \right)}}\end{matrix}\right]\end{split}\]

But it is best to use the LUsolve() method to perform an LU decomposition Gaussian-Elimination to solve the system, especially as the dimension of \(\mathbf{A}\) grows:

A.LUsolve(b)

\[\begin{split}\displaystyle \left[\begin{matrix}\frac{- \frac{\omega \log{\left(t \right)}}{\log{\left(f{\left(t \right)} \right)}} - 100 - \frac{- \frac{\omega^{2} \left(- \frac{\omega \log{\left(t \right)}}{\log{\left(f{\left(t \right)} \right)}} - 100\right)}{\sin{\left(f{\left(t \right)} \right)} \cos{\left(2 f{\left(t \right)} \right)}} - \omega - f^{3}{\left(t \right)}}{- \frac{\omega^{2}}{\sin{\left(f{\left(t \right)} \right)} \cos{\left(2 f{\left(t \right)} \right)}} + f{\left(t \right)}}}{\sin{\left(f{\left(t \right)} \right)} \cos{\left(2 f{\left(t \right)} \right)}}\\\frac{- \frac{\omega^{2} \left(- \frac{\omega \log{\left(t \right)}}{\log{\left(f{\left(t \right)} \right)}} - 100\right)}{\sin{\left(f{\left(t \right)} \right)} \cos{\left(2 f{\left(t \right)} \right)}} - \omega - f^{3}{\left(t \right)}}{- \frac{\omega^{2}}{\sin{\left(f{\left(t \right)} \right)} \cos{\left(2 f{\left(t \right)} \right)}} + f{\left(t \right)}}\end{matrix}\right]\end{split}\]

Warning

This method of solving symbolic linear systems is fast, but it can give incorrect answers for:

1. expressions that are not acutally linear in the variables the Jacobian is taken with respect to

2. \(\mathbf{A}\) matrix entries that would evaluate to zero if simplified or specific numerical values are provided

So only use this method if you are sure your equations are linear and if your \(\mathbf{A}\) matrix is made up of complex expressions, watch out for nan results after lambdifying. solve() and linsolve() can also solve linear systems and they check for linearity and properties of the A matrix. The cost is that they can be extremely slow for large expressions (which we will have in this book).

Exercise

Solve the following equations for all of the \(L\)’s and then use lambdify() to evaluate the solution for \(F_1=13\) and \(F_2=32\).

(10)\[\begin{split}-L_1 + L_2 - L_3/\sqrt{2} = & 0 \\ L_3/\sqrt{2} + L_4 = & F_1 \\ -L_2 - L_5/\sqrt{2} = & 0 \\ L_5/\sqrt{2} = & F_2 \\ L_5/\sqrt{2} + L_6 = & 0 \\ -L_4 -L_5/\sqrt{2} = & 0\end{split}\]

Solution

L1, L2, L3, L4, L5, L6, F1, F2 = sm.symbols('L1, L2, L3, L4, L5, L6, F1, F2')

exprs = sm.Matrix([
    -L1 + L2 - L3/sm.sqrt(2),
    L3/sm.sqrt(2) + L4 - F1,
    -L2 - L5/sm.sqrt(2),
    L5/sm.sqrt(2) - F2,
    L5/sm.sqrt(2) + L6,
    -L4 -L5/sm.sqrt(2),
])
exprs

\[\begin{split}\displaystyle \left[\begin{matrix}- L_{1} + L_{2} - \frac{\sqrt{2} L_{3}}{2}\\- F_{1} + \frac{\sqrt{2} L_{3}}{2} + L_{4}\\- L_{2} - \frac{\sqrt{2} L_{5}}{2}\\- F_{2} + \frac{\sqrt{2} L_{5}}{2}\\\frac{\sqrt{2} L_{5}}{2} + L_{6}\\- L_{4} - \frac{\sqrt{2} L_{5}}{2}\end{matrix}\right]\end{split}\]

unknowns = sm.Matrix([L1, L2, L3, L4, L5, L6])

coef_mat = exprs.jacobian(unknowns)
rhs = -exprs.xreplace(dict(zip(unknowns, [0]*6)))

sol = coef_mat.LUsolve(rhs)

sm.Eq(unknowns, sol)

\[\begin{split}\displaystyle \left[\begin{matrix}L_{1}\\L_{2}\\L_{3}\\L_{4}\\L_{5}\\L_{6}\end{matrix}\right] = \left[\begin{matrix}- F_{1} - 2 F_{2}\\- F_{2}\\\sqrt{2} \left(F_{1} + F_{2}\right)\\- F_{2}\\\sqrt{2} F_{2}\\- F_{2}\end{matrix}\right]\end{split}\]

eval_sol = sm.lambdify((F1, F2), sol)
eval_sol(13, 32)

array([[-77.        ],
       [-32.        ],
       [ 63.63961031],
       [-32.        ],
       [ 45.254834  ],
       [-32.        ]])

Simplification

The above result from LUsolve() is a bit complicated. Reproduced here:

a1, a2 = sm.symbols('a1, a2')
exprs = sm.Matrix([
    [a1*sm.sin(f(t))*sm.cos(2*f(t)) + a2 + omega/sm.log(f(t), t) + 100],
    [a1*omega**2 + f(t)*a2 + omega + f(t)**3],
])
A = exprs.jacobian([a1, a2])
b = -exprs.xreplace({a1: 0, a2: 0})
sol = A.LUsolve(b)

SymPy has some functionality for automatically simplifying symbolic expressions. The function simplify() will attempt to find a simpler version:

sm.simplify(sol)

\[\begin{split}\displaystyle \left[\begin{matrix}\frac{- \omega f{\left(t \right)} \log{\left(t \right)} + \omega \log{\left(f{\left(t \right)} \right)} + f^{3}{\left(t \right)} \log{\left(f{\left(t \right)} \right)} - 100 f{\left(t \right)} \log{\left(f{\left(t \right)} \right)}}{\left(- \omega^{2} + f{\left(t \right)} \sin{\left(f{\left(t \right)} \right)} \cos{\left(2 f{\left(t \right)} \right)}\right) \log{\left(f{\left(t \right)} \right)}}\\\frac{- \omega^{2} \left(\omega \log{\left(t \right)} + 100 \log{\left(f{\left(t \right)} \right)}\right) + \left(\omega + f^{3}{\left(t \right)}\right) \log{\left(f{\left(t \right)} \right)} \sin{\left(f{\left(t \right)} \right)} \cos{\left(2 f{\left(t \right)} \right)}}{\left(\omega^{2} - f{\left(t \right)} \sin{\left(f{\left(t \right)} \right)} \cos{\left(2 f{\left(t \right)} \right)}\right) \log{\left(f{\left(t \right)} \right)}}\end{matrix}\right]\end{split}\]

But you’ll have the best luck at simplifying if you use simplification functions that target the type of expression you have. The trigsimp() function only attempts trigonometric simplifications, for example:

trig_expr = sm.cos(omega)**2 + sm.sin(omega)**2
trig_expr

\[\displaystyle \sin^{2}{\left(\omega \right)} + \cos^{2}{\left(\omega \right)}\]

sm.trigsimp(trig_expr)

\[\displaystyle 1\]

Warning

Only attempt simplification on expressions that are several lines of text. Larger expressions become increasingly computationally intensive to simplify and there is generally no need to do so.

As mentioned earlier, SymPy represents expressions as trees. Symbolic expressions can also be represented as directed acyclic graphs that contain only one node for each unique expression (unlike SymPy’s trees which may have the same expression in more than one node). These unique expressions, or “common subexpressions”, can be found with the cse() function. This function will provide a simpler form of the equations that minimizes the number of operations to compute the answer. We can count the number of basic operations (additions, multiplies, etc.) using count_ops():

sm.count_ops(sol)

\[\displaystyle 79\]

We can simplify with cse():

substitutions, simplified = sm.cse(sol)

The substitutions variable contains a list of tuples, where each tuple has a new intermediate variable and the sub-expression it is equal to.

substitutions[0]

\[\displaystyle \left( x_{0}, \ f{\left(t \right)}\right)\]

The Eq() class with tuple unpacking (*) can be used to display these tuples as equations:

sm.Eq(*substitutions[0])

\[\displaystyle x_{0} = f{\left(t \right)}\]

sm.Eq(*substitutions[1])

\[\displaystyle x_{1} = \frac{1}{\sin{\left(x_{0} \right)} \cos{\left(2 x_{0} \right)}}\]

sm.Eq(*substitutions[2])

\[\displaystyle x_{2} = \omega^{2} x_{1}\]

sm.Eq(*substitutions[4])

\[\displaystyle x_{4} = \frac{- \omega - x_{0}^{3} + x_{2} x_{3}}{x_{0} - x_{2}}\]

The simplified variable contains the simplified expression, made up of the intermediate variables.

simplified[0]

\[\begin{split}\displaystyle \left[\begin{matrix}x_{1} \left(- x_{3} - x_{4}\right)\\x_{4}\end{matrix}\right]\end{split}\]

We can count the number of operations of the simplified version:

num_ops = sm.count_ops(simplified[0])
for sub in substitutions:
    num_ops += sm.count_ops(sub[1])
num_ops

\[\displaystyle 22\]

Exercise

lambdify() has an optional argument cse=True|False that applies common subexpression elimination internally to simplify the number of operations. Differentiate the base_expr with respect to x 10 times to generate a very long expression. Create two functions using lambdify(), one with cse=True and one with cse=False. Compare how long it takes to numerically evaluate the resulting functions using the %timeit magic.

a, b, c, x, y, z = sm.symbols('a, b, c, x, y, z')
base_expr = a*sm.sin(x*x + b*sm.cos(x*y) + c*sm.sin(x*z))

Solution

Differentiate 10 times:

long_expr = base_expr.diff(x, 10)

Create the numerical functions:

eval_long_expr = sm.lambdify((a, b, c, x, y, z), long_expr)
eval_long_expr_cse = sm.lambdify((a, b, c, x, y, z), long_expr, cse=True)

Now time each function:

%%timeit
eval_long_expr(1.0, 2.0, 3.0, 4.0, 5.0, 6.0)

261 μs ± 1.13 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

%%timeit
eval_long_expr_cse(1.0, 2.0, 3.0, 4.0, 5.0, 6.0)

16.5 μs ± 31.8 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

Learn more

This section only scratches the surface of what SymPy can do. The presented concepts are the basic ones needed for this book, but getting more familiar with SymPy and what it can do will help. I recommend doing the SymPy Tutorial. The “Gotchas” section is particularly helpful for common mistakes when using SymPy. The tutorial is part of the SymPy documentation https://docs.sympy.org, where you will find general information on SymPy.

The tutorial is also available on video:

If you want to ask a question about using SymPy (or search to see if someone else has asked your question), you can do so at the following places:

• SymPy mailing list: Ask questions via email.

• SymPy Github Discussions: Ask questions via Github.

• Stackoverflow: Ask and search questions on the most popular coding Q&A website.