Symbolic Math in MATLAB with SymPy
Part 10 of the Python Is The Ultimate MATLAB Toolbox series.
Introduction
SymPy is a Python-based computer algebra system with capabilities similar to Maple, Maxima, MATLAB’s Symbolic Math Toolbox, SageMath, and to a lesser extent, Mathematica. The most natural way to work with mathematical expressions symbolically in MATLAB is with the Symbolic Math Toolbox. If you don’t have access to this toolbox and you’re willing to include Python in your MATLAB workflow, SymPy may be an adequate alternative. SymPy is included with the recommended Anaconda Python distribution.
SymPy’s capabilities are vast—the PDF version of its recent v1.11.1 documentation is a staggering 3,184 pages—and this post merely covers the basics: defining symbolic variables, integrating, taking derivatives, solving simultaneous equations, summing series, generating LaTeX markup, and generating source code in MATLAB and other languages.
To give you a sense of SymPy’s depth, here are a few of its capabilities that aren’t covered in this article: linear, time-invariant (LTI) controls, combinatorics, quaternions, sets, tensor products and contraction, Pauli algebra, hydrogen wavefunctions, quantum harmonic oscillators, gamma matrices, Gaussian optics, and continuum mechanics.
- Defining symbolic variables and functions
- Solving equations
- Solution as LaTeX markup
- Solution as MATLAB / Python / C / Fortran / etc source code
- Definite and indefinite integrals
- Derivatives
- Summing series
- Laplace transforms
Defining symbolic variables and functions
The first steps to using SymPy are importing its various modules
and defining symbolic variables and functions.
Symbolic variables can be defined three ways: 1) imported
from sympy.abc
, 2) created with a call to sympy.symbols()
,
or 3) created with a call to sympy.Symbol()
.
The module sympy.abc
has a collection of predefined symbolic
variables that include all lower- and uppercase
letters and words defining Greek letters such as omega
.
Multi-character variables can be defined with
the sympy.symbols()
or sympy.Symbol()
functions. The latter
allows one to add properties to a variable indicating, for
example, that it is an integer, or
its value is always positive, or that it is a constant,
or must be real.
The code below defines symbolic variables x
, y
, r
, phi
,
u1
, v1
, and Length
:
Python:
from sympy import symbols, Symbol
from sympy.abc import x, y, r, phi # method 1
u1, v1 = symbols('u1 v1') # method 2
Length = Symbol('Length', positive=True, constant=True, real=True) # method 3
The constants $\pi$, $i$, $e$, and $\infty$ must be imported to be used
symbolically. $i$ and $e$ exist as uppercase I
and E
, and
$\infty$ is oo
:
Python:
In : from sympy import Float, pi, I, E, oo
In : pi, I, E, oo
Out: (pi, I, E, oo)
In : Float(pi), Float(E), Float(oo)
Out: (3.14159265358979, 2.71828182845905, oo)
Float(I)
doesn’t appear above because yields a TypeError
.
Unlike $\pi$ and $e$, infinity is not expanded to its IEEE 754 floating
point representation even though NumPy has an explicit term for it: np.inf
.
Symbolic functions
In addition to symbolic variables, functions that operate on symbolic variables must also be explicitly requested from SymPy. Symbolic-capable functions include forward and inverse trignometric, forward and inverse hyperbolic, logarithmic, and many others.
For example if your work involves symbolic computations with sine, cosine, and natural logarithm functions, load the symbolic versions of these functions from SymPy with:
Python:
from sympy import sin, cos, ln
Examples of equivalent MATLAB imports appear in the integration and derivation sections below.
Solving equations
We’ll use the symbols defined above in some equations and have SymPy solve them. Here’s an easy one to start: given the equation for a circle,
then supplying the left hand side to
solve()
, along with the variable we want to solve for, namely x
:
Python:
In : from sympy.abc import x, y, r
In : from sympy import solve
In : solve(x**2 + y**2 - r**2, x)
Out: [-sqrt((r - y)*(r + y)), sqrt((r - y)*(r + y))]
Python:
In : from sympy.abc import x, y, r
In : from sympy import solve, expand
In : sol_minus, sol_plus = solve(x**2 + y**2 - r**2, x)
In : expand(sol_minus), expand(sol_plus)
Out: (-sqrt(r**2 - y**2), sqrt(r**2 - y**2))
… which raises the question of what does it mean for an equation to be “simplified” or “expanded”? SymPy’s opinion in this case is the opposite of my understanding of these terms.
The second way to solve an equation is to create an
explicit equation object
with the Eq()
function then pass the equation to solve()
.
Eq()
takes two input arguments which
are the left and right sides of the equation:
Python:
In : from sympy.abc import x, y, r
In : from sympy import Eq
In : circle_eqn = Eq(x**2 + y**2, r**2) # x^2 + y^2 = r^2
In : solve(circle_eqn, x)
Out: [-sqrt((r - y)*(r + y)), sqrt((r - y)*(r + y))]
Intersection of a circle and line in Python
So far we haven’t really solved anything, we’ve just reordered terms of an equation. Next we’ll add the equation of a line to the mix and have SymPy figure out the points of intersection. The equations are
This is a nonlinear system of equations so we’ll need SymPy’s
nonlinsolve()
function:
Python:
In : from sympy import Eq, nonlinsolve
In : from sympy.abc import x, y, r, m, b
In : circle_eqn = Eq(x**2 + y**2, r**2) # x^2 + y^2 = r^2
In : line_eqn = Eq(y, m*x + b) # y = m x + b
In : nonlinsolve([circle_eqn, line_eqn], [x, y])
Out: {((-b + b/(m**2 + 1) - m*sqrt(-b**2 + m**2*r**2 + r**2)/(m**2 + 1))/m,
b/(m**2 + 1) - m*sqrt(-b**2 + m**2*r**2 + r**2)/(m**2 + 1)),
((-b + b/(m**2 + 1) + m*sqrt(-b**2 + m**2*r**2 + r**2)/(m**2 + 1))/m,
b/(m**2 + 1) + m*sqrt(-b**2 + m**2*r**2 + r**2)/(m**2 + 1))}
SymPy returns two pairs of $(x,y)$ values representing the points of intersection—if these exist, of course. The actual values of $b$, $m$, and $r$ will determine whether the solution has real or imaginary values.
In Solve Systems of Nonlinear Equations with scipy.optimize.root()
I show how SciPy can find numeric solutions to
systems of nonlinear equations.
How does SymPy fare on the intersection of four spheres problem
described there?
Unfortunately, not well. The program below assembles the problem
but nonlinearsolve()
doesn’t return after running for more than an hour.
Python:
|
|
Intersection of a circle and line in MATLAB
The MATLAB version of the circle/line intersection problem looks like this:
MATLAB (2020b and newer):
|
|
Note: In MATLAB, exponents for squares, cubes, and so on should be explicitly cast to integers in the SymPy expression. If this is not done, the notation “
^2
” given to SymPy will use the default MATLAB type ofdouble
which greatly complicates — and often defeats — the SymPy evaluator.
The solution in MATLAB appears with the full complement of attributes in addition to just the equations we’re after. This is a truncated portion of the output:
MATLAB:
>> soln
Python FiniteSet with properties:
args: [1x2 py.tuple]
assumptions0: [1x1 py.dict]
boundary: [1x1 py.sympy.sets.sets.FiniteSet]
:
is_transcendental: [1x1 py.NoneType]
is_zero: [1x1 py.NoneType]
measure: [1x1 py.int]
{((-b + b/(m**2 + 1) - m*sqrt(-b**2 + m**2*r**2 + r**2)/(m**2 + 1))/m,
b/(m**2 + 1) - m*sqrt(-b**2 + m**2*r**2 + r**2)/(m**2 + 1)),
((-b + b/(m**2 + 1) + m*sqrt(-b**2 + m**2*r**2 + r**2)/(m**2 + 1))/m,
b/(m**2 + 1) + m*sqrt(-b**2 + m**2*r**2 + r**2)/(m**2 + 1))}
The string containing just the two $(x,y)$ pairs is in soln.char
:
MATLAB:
>> soln.char
'{((-b + b/(m**2 + 1) - m*sqrt(-b**2 + m**2*r**2 + r**2)/(m**2 + 1))/m,
b/(m**2 + 1) - m*sqrt(-b**2 + m**2*r**2 + r**2)/(m**2 + 1)),
((-b + b/(m**2 + 1) + m*sqrt(-b**2 + m**2*r**2 + r**2)/(m**2 + 1))/m,
b/(m**2 + 1) + m*sqrt(-b**2 + m**2*r**2 + r**2)/(m**2 + 1))}'
Solution as LaTeX
The SymPy text results we’ve seen so far look awkward. If you want to publish the results they’ll look much more elegant rendered with LaTeX. Simple expressions are easy to write manually but SymPy’s output can get lengthy so its LaTeX generator comes in handy:
Python:
In : from sympy.printing.latex import print_latex
In : from sympy import Eq, nonlinsolve
In : from sympy.abc import x, y, r, m, b
In : circle_eqn = Eq(x**2 + y**2, r**2) # x^2 + y^2 = r^2
In : line_eqn = Eq(y, m*x + b) # y = m x + b
In : soln = nonlinsolve([circle_eqn, line_eqn], [x, y])
In : print_latex(soln)
The last line produces this output
\left\{\left( \frac{- b + \frac{b}{m^{2} + 1} - \frac{m \sqrt{- b^{2} + m^{2} r^{2} + r^{2}}}{m^{2} + 1}}{m}, \ \frac{b}{m^{2} + 1} - \frac{m \sqrt{- b^{2} + m^{2} r^{2} + r^{2}}}{m^{2} + 1}\right), \left( \frac{- b + \frac{b}{m^{2} + 1} + \frac{m \sqrt{- b^{2} + m^{2} r^{2} + r^{2}}} {m^{2} + 1 }}{m}, \ \frac{b}{m^{2} + 1} + \frac{m \sqrt{- b^{2} + m^{2} r^{2} + r^{2}}}{m^{2} + 1}\right)\right\}
which renders as
Solution as MATLAB / Python / C / Fortran / etc source code
In addition to LaTeX, SymPy can express solutions as computer
source code in a dozen languages. The
sympy.printing
module has these emitters:
Language | sympy.printing function |
---|---|
C | ccode() |
C++ | cxxcode() |
Fortran | fcode() |
GLSL | glsl_code() |
JavaScript | jscode() |
Julia | julia_code() |
Maple | maple_code() |
Mathematica | mathematica_code() |
MATLAB | octave_code() |
Python | pycode() |
R | rcode() |
Rust | rust_code() |
The code blocks below repeat the circle/line intersection
this time emitting MATLAB source code instead of LaTeX.
In MATLAB, the call to octave_code()
returns a Python string
so this is converted to a MATLAB string before printing.
Python:
|
|
MATLAB (2020b and newer):
|
|
Both versions produce this output:
"{{(-b + b./(m.^2 + 1) - m.*sqrt(-b.^2 + m.^2.*r.^2 + r.^2)./(m.^2 + 1))./m,
b./(m.^2 + 1) - m.*sqrt(-b.^2 + m.^2.*r.^2 + r.^2)./(m.^2 + 1)},
{(-b + b./(m.^2 + 1) + m.*sqrt(-b.^2 + m.^2.*r.^2 + r.^2)./(m.^2 + 1))./m,
b./(m.^2 + 1) + m.*sqrt(-b.^2 + m.^2.*r.^2 + r.^2)./(m.^2 + 1)}}"
Replacing octave_code()
with jscode()
produces much lengthier
JavaScript code:
{((-b + b/(Math.pow(m, 2) + 1) - m*Math.sqrt(-Math.pow(b, 2) +
Math.pow(m, 2)*Math.pow(r, 2) + Math.pow(r, 2))/(Math.pow(m, 2) + 1))/m,
b/(Math.pow(m, 2) + 1) - m*Math.sqrt(-Math.pow(b, 2) +
Math.pow(m, 2)*Math.pow(r, 2) + Math.pow(r, 2))/(Math.pow(m, 2) + 1)),
((-b + b/(Math.pow(m, 2) + 1) + m*Math.sqrt(-Math.pow(b, 2) +
Math.pow(m, 2)*Math.pow(r, 2) + Math.pow(r, 2))/(Math.pow(m, 2) + 1))/m,
b/(Math.pow(m, 2) + 1) + m*Math.sqrt(-Math.pow(b, 2) +
Math.pow(m, 2)*Math.pow(r, 2) + Math.pow(r, 2))/(Math.pow(m, 2) + 1))}
Definite and indefinite integrals
Integration is done with SymPy’s integrate()
function. Here’s
a tricky one:
In addition to the integrate()
function we’ll also need to
import symbolic versions of the sine and cosine functions.
Constraining $n$ and $x$ to be real, and $n$ to be positive
yields a cleaner result.
The Python and MATLAB solutions are
Python:
|
|
MATLAB:
|
|
Note: In the MATLAB solution, SymPy versions of $\pi$ and sine and cosine functions are imported with a leading uppercase character to make them distinct from MATLAB’s
pi
,sin()
, andcos()
. If you importpi
,sin
, andcos
from SymPy using the identical names they will be SymPy objects that cannot be used in conventional MATLAB expressions.
The text version of the solution, minus the obligatory constant, from Python is
7*sin(pi*n*x)*sin(3*pi*n*x/2)**2/(16*pi*n) +
9*sin(pi*n*x)*cos(3*pi*n*x/2)**2/(16*pi*n) + sin(pi*n*x)/(pi*n) -
3*sin(3*pi*n*x/2)*cos(pi*n*x)*cos(3*pi*n*x/2)/(8*pi*n)
which renders like this in LaTeX:
The MATLAB result looks quite different
0.4375*sin(pi*n*x)*sin(3*pi*n*x/2)**2/(pi*n) +
0.5625*sin(pi*n*x)*cos(3*pi*n*x/2)**2/(pi*n) + 1.0*sin(pi*n*x)/(pi*n) -
0.375*sin(3*pi*n*x/2)*cos(pi*n*x)*cos(3*pi*n*x/2)/(pi*n)
but, as shown with numeric evaluation, both yield the same results.
Numeric evaluation of symbolic results
One way to convince ourselves that two starkly different equations
are equivalent is to substitute numeric values in for the variables.
SymPy equations have a .evalf()
method to facilitate this
substitution. The method accepts a dictionary of substitutions
and returns a numeric value.
The last two lines of the Python code above,
Python:
at_x_n = soln.evalf(subs={x : 1.5, n : 1})
print(f'solution at x=1.5, n=1 is {at_x_n}')
solution at x=1.5, n=1 is -0.477464829275686
I’ve not been able to replicate .evalf()'s
subs=
expression in MATLAB, however.
Instead, to compare the results of the two integrals,
I generated MATLAB source code for each and substituted
arrays of values and confirmed both Python and MATLAB
solutions agree.
Definite integrals
Definite integrals are defined by supplying lower and upper bounds after the variable of integration. The value of the same integral shown above going from $x = 0$ to $x = \frac{1}{\pi}$ for $n = 1$, in other words
is solved like this in Python and MATLAB:
Python:
|
|
MATLAB:
|
|
Once again, the Python and MATLAB symbolic solutions differ but both evaluate to the same number:
Python:
MATLAB:
In this case .evalf()
works in MATLAB because there are no
variables to substitute.
Derivatives
Derivatives are computed with the diff()
function. If we start
with a function
we can get $\frac{dy}{dx}$ with
Python:
|
|
MATLAB:
|
|
both of which give this LaTeX representation of the solution:
Partial derivatives
diff()
can also return partial derivatives. Say you have a
complicated surface such as
and you want to compute the matrix of second order derivatives
you could do it with
Python:
|
|
MATLAB:
|
|
which return different arrangements of the same result:
Python:
MATLAB:
Summing series
SymPy has support for Fourier series, power series, series summation, and series expansion. This example shows how the closed-form summation
can be found:
Python:
|
|
MATLAB:
|
|
which give an answer of
Laplace transforms
The laplace_transform()
function in SymPy resembles
the laplace()
function in MATLAB’s Symbolic Toolbox.
Here we’ll use it to compute the transform of
Python:
|
|
MATLAB:
|
|
which produces LaTeX that renders as
Table of Contents | Previous: Solve simultaneous nonlinear equations Next : Accelerate MATLAB with Python and Numba