Recent Posts (page 3 / 8)

by Al Danial

Symbolic Math in MATLAB with SymPy

Python MATLAB CAS symbolic math SymPy integrate Laplace transform sum series derivatives integrals

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

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,

$$ x^2 + y^2 = r^2 $$
express the circle in terms of $x$. There are two ways to set up the problem for SymPy. The first involves recasting the equation to equal zero, as in
$$ x^2 + y^2 - r^2 = 0 $$

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))]
While the answer is correct, it doesn’t resemble the expected form of

$$ x = \pm \sqrt{ r^2 - y^2 } $$
We can get the expected form by explicitly expanding the two parts of the solution:

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

$$ \begin{align*} x^2 + y^2 &= r^2 \\ y &= m x + b \end{align*} $$

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: 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

#!/usr/bin/env python
# file: four_spheres.py

# Don't try this at home; the nonlinsolve() step hangs.

from sympy import Eq, nonlinsolve, symbols
from sympy.abc import x, y, z, c, e
A1, A2, A3, A4 = symbols('A1 A2 A3 A4')
B1, B2, B3, B4 = symbols('B1 B2 B3 B4')
C1, C2, C3, C4 = symbols('C1 C2 C3 C4')
Dt1, Dt2, Dt3, Dt4 = symbols('Dt1 Dt2 Dt3 Dt4')
s1 = Eq( (x-A1)**2 + (y-B1)**2 + (z-C1)**2 , (c*(Dt1 + e) )**2 )
s2 = Eq( (x-A2)**2 + (y-B2)**2 + (z-C2)**2 , (c*(Dt2 + e) )**2 )
s3 = Eq( (x-A3)**2 + (y-B3)**2 + (z-C3)**2 , (c*(Dt3 + e) )**2 )
s4 = Eq( (x-A4)**2 + (y-B4)**2 + (z-C4)**2 , (c*(Dt4 + e) )**2 )
soln = nonlinsolve([s1, s2, s3, s4], [x, y, z, e])  # hangs

Intersection of a circle and line in MATLAB

The MATLAB version of the circle/line intersection problem looks like this:

MATLAB (2020b and newer): 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

sympy = py.importlib.import_module('sympy');
abc = py.importlib.import_module('sympy.abc');
x = abc.x;
y = abc.y;
r = abc.r;
m = abc.m;
b = abc.b;
Two = int64(2);
circle_eqn = sympy.Eq(x^Two + y^Two, r^Two); % x^2 + y^2 = r^2
line_eqn   = sympy.Eq(y, m*x + b);           % y = m*x + b
soln = sympy.nonlinsolve(py.list({circle_eqn, line_eqn}), ...
                         py.list({x, y}))

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 of double 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

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

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: 

1
2
3
4
5
6
7
8
9

#!/usr/bin/env python
from sympy import Eq, nonlinsolve
from sympy.abc import x, y, r, m, b
from sympy.printing import octave_code
circle_eqn = Eq(x**2 + y**2, r**2) # x^2 + y^2 = r^2
line_eqn = Eq(y, m*x + b)          # y = m x + b
soln = nonlinsolve([circle_eqn, line_eqn], [x, y])
print(soln)
print(octave_code(soln))

MATLAB (2020b and newer): 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

sympy = py.importlib.import_module('sympy');
abc = py.importlib.import_module('sympy.abc');
sym_code = py.importlib.import_module('sympy.printing');
x = abc.x;
y = abc.y;
r = abc.r;
m = abc.m;
b = abc.b;
Two = int64(2);
circle_eqn = sympy.Eq(x^Two + y^Two, r^Two); % x^2 + y^2 = r^2
line_eqn   = sympy.Eq(y, m*x + b);           % y = m*x + b
soln = sympy.nonlinsolve(py.list({circle_eqn, line_eqn}), ...
                         py.list({x, y}));
fprintf('%s\n', string(sym_code.octave_code(soln)) )

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:

$$ \int \cos(\pi n x) \left( 1 + \sin^2 \left(\frac{3 \pi n x}{2} \right) \right) dx $$

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: 

1
2
3
4
5
6
7
8
9

#!/usr/bin/env python
from sympy import pi, sin, cos, integrate, Symbol
from sympy.printing.latex import print_latex
from sympy.printing import octave_code
n = Symbol('n', real=True, positive=True, integer=True)
x = Symbol('x', real=True)
soln = integrate(cos(pi*n*x)*(1 + sin(3*pi*n*x/2)**2), (x))
at_x_n = soln.evalf(subs={x : 1.5, n : 1})
print(f'solution at x=1.5, n=1 is {at_x_n}')

MATLAB: 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

sympy = py.importlib.import_module('sympy');
printing = py.importlib.import_module('sympy.printing');
n = sympy.Symbol('n', pyargs('real',py.True, ...
                             'positive',py.True,'integer',py.True));
x = sympy.Symbol('x', pyargs('real',py.True));
Pi = sympy.pi;   % uppercase to be distinct from matlab's pi
Sin = sympy.sin; % uppercase to be distinct from matlab's sin()
Cos = sympy.cos; % uppercase to be distinct from matlab's cos()
soln = sympy.integrate(Cos(Pi*n*x)*(1 + Sin(3*Pi*n*x/2)^int64(2)), ...
                       py.tuple({x}));
fprintf('%s\n', soln.char)
printing.print_latex(soln)

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(), and cos(). If you import pi, sin, and cos 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:

$$ \frac{7 \sin{\left(\pi n x \right)} \sin^{2}{\left(\frac{3 \pi n x}{2} \right)}}{16 \pi n} + \frac{9 \sin{\left(\pi n x \right)} \cos^{2}{\left(\frac{3 \pi n x}{2} \right)}}{16 \pi n} + \frac{\sin{\left(\pi n x \right)}}{\pi n} - \frac{3 \sin{\left(\frac{3 \pi n x}{2} \right)} \cos{\left(\pi n x \right)} \cos{\left(\frac{3 \pi n x}{2} \right)}}{8 \pi n} $$

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)
$$ \frac{0.4375 \sin{\left(\pi n x \right)} \sin^{2}{\left(\frac{3 \pi n x}{2} \right)}}{\pi n} + \frac{0.5625 \sin{\left(\pi n x \right)} \cos^{2}{\left(\frac{3 \pi n x}{2} \right)}}{\pi n} + \frac{1.0 \sin{\left(\pi n x \right)}}{\pi n} - \\ \frac{0.375 \sin{\left(\frac{3 \pi n x}{2} \right)} \cos{\left(\pi n x \right)} \cos{\left(\frac{3 \pi n x}{2} \right)}}{\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}')

print

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

$$ \int_0^{\frac{1}{\pi}} \cos(\pi x) \left( 1 + \sin^2 \left(\frac{3 \pi x}{2} \right) \right) dx $$

is solved like this in Python and MATLAB:

Python: 

1
2
3
4
5
6
7

#!/usr/bin/env python
from sympy import pi, sin, cos, integrate, Symbol, simplify
from sympy.printing.latex import print_latex
from sympy.printing import octave_code
x = Symbol('x', real=True)
soln = integrate(cos(pi*x)*(1 + sin(3*pi*x/2)**2), (x, 0, 1/pi))
print_latex(simplify(soln))

MATLAB: 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

sympy = py.importlib.import_module('sympy');
printing = py.importlib.import_module('sympy.printing');
x = sympy.Symbol('x', pyargs('real',py.True));
Pi = sympy.pi;   % uppercase to be distinct from matlab's pi
Sin = sympy.sin; % uppercase to be distinct from matlab's sin()
Cos = sympy.cos; % uppercase to be distinct from matlab's cos()
soln = sympy.integrate(Cos(Pi*x)*(1 + Sin(3*Pi*x/2)^int64(2)), ...
                       py.tuple({x, 0, 1/Pi}));
printing.print_latex(soln)
fprintf('%s\n', soln.evalf().char)

Once again, the Python and MATLAB symbolic solutions differ but both evaluate to the same number:

Python:

$$ \frac{- 2 \sin{\left(2 \right)} - \sin{\left(4 \right)} + 24 \sin{\left(1 \right)}}{16 \pi} = 0.380649112305801 $$

MATLAB:

$$ \frac{1.19584445481538}{\pi} = 0.380649112305801 $$

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

$$ y = \sin^2(\pi x) e^{- i x} $$

we can get $\frac{dy}{dx}$ with

Python: 

1
2
3
4
5
6

#!/usr/bin/env python
from sympy import I, pi, sin, exp, diff
from sympy.abc import x, y
from sympy.printing.latex import print_latex
soln = diff(sin(pi*x)**2 * exp(-I*x), x)
print_latex(soln)

MATLAB: 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

sympy = py.importlib.import_module('sympy');
abc = py.importlib.import_module('sympy.abc');
printing = py.importlib.import_module('sympy.printing');
x = abc.x;
y = abc.y;
I = sympy.I;
Pi = sympy.pi;   % uppercase to be distinct from matlab's pi
Sin = sympy.sin; % uppercase to be distinct from matlab's sin()
Exp = sympy.exp; % uppercase to be distinct from matlab's exp()
Two = int64(2);
soln = sympy.diff(Sin(Pi*x)^Two * Exp(-I*x), x);
fprintf('%s\n', soln.char)
printing.print_latex(soln)

both of which give this LaTeX representation of the solution:

$$ \frac{dy}{dx} = - i e^{- i x} \sin^{2}{\left(\pi x \right)} + 2 \pi e^{- i x} \sin{\left(\pi x \right)} \cos{\left(\pi x \right)} $$

Partial derivatives

diff() can also return partial derivatives. Say you have a complicated surface such as

$$ z(x,y) = 4 x^4 - 17 x^3 y + 3 x^2 y^2 - y^3 - y^4 $$

and you want to compute the matrix of second order derivatives

$$ \begin{bmatrix} \frac{\partial^2 z}{\partial^2 x} & \frac{\partial^2 z}{\partial x \partial y} \\ \frac{\partial^2 z}{\partial y \partial x} & \frac{\partial^2 z}{\partial^2 y} \end{bmatrix} $$

you could do it with

Python: 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

#!/usr/bin/env python
from sympy import diff
from sympy.abc import x, y, z
from sympy.printing.latex import print_latex
z = 4*x**4 - 17*x**3*y + 3*x**2*y**2 - y**3 - y**4
dzz_dxx = diff(z, x, x)
dzz_dxy = diff(z, x, y)
dzz_dyx = diff(z, y, x)
dzz_dyy = diff(z, y, y)
print_latex(dzz_dxx)
print_latex(dzz_dxy)
print_latex(dzz_dyx)
print_latex(dzz_dyy)

MATLAB: 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

sympy = py.importlib.import_module('sympy');
abc = py.importlib.import_module('sympy.abc');
printing = py.importlib.import_module('sympy.printing');
x = abc.x;
y = abc.y;
z = abc.z;
E2 = int64(2);
E3 = int64(3);
E4 = int64(4);
z = 4*x^E4 - 17*x^E3*y + 3*x^E2*y^E2 - y^E3 - y^E4;
dzz_dxx = diff(z, x, x);
dzz_dxy = diff(z, x, y);
dzz_dyx = diff(z, y, x);
dzz_dyy = diff(z, y, y);
printing.print_latex(dzz_dxx)
printing.print_latex(dzz_dxy)
printing.print_latex(dzz_dyx)
printing.print_latex(dzz_dyy)

which return different arrangements of the same result:

Python:

$$ \begin{bmatrix} 6 \cdot \left(8 x^{2} - 17 x y + y^{2}\right) & 3 x \left(- 17 x + 4 y\right) \\ 3 x \left(- 17 x + 4 y\right) & 6 \left(x^{2} - 2 y^{2} - y\right) \end{bmatrix} $$

MATLAB:

$$ \begin{bmatrix} 48.0 x^{2} - 102.0 x y + 6.0 y^{2} & x \left(- 51.0 x + 12.0 y\right) \\ x \left(- 51.0 x + 12.0 y\right) & 6.0 x^{2} - 12 y^{2} - 6 y \end{bmatrix} $$

Summing series

SymPy has support for Fourier series, power series, series summation, and series expansion. This example shows how the closed-form summation

$$ \displaystyle\sum_{k=1}^{N+1} \frac{e^{ik\pi}}{N^2} $$

can be found:

Python: 

1
2
3
4
5
6
7

#!/usr/bin/env python
from sympy import I, Sum, exp, pi
from sympy.abc import k, N
from sympy.printing.latex import print_latex
s = Sum( exp(I*k*pi)/N**2, (k, 1, N+1))
soln = s.doit()
print_latex(soln)

MATLAB: 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13

sympy = py.importlib.import_module('sympy');
abc = py.importlib.import_module('sympy.abc');
printing = py.importlib.import_module('sympy.printing');
k = abc.k;
N = abc.N;
I = sympy.I;
Pi  = sympy.pi;  % uppercase to be distinct from matlab's pi
Exp = sympy.exp; % uppercase to be distinct from matlab's exp()
Sum = sympy.Sum;
Two = int64(2);
s = Sum( Exp(I*k*Pi)/(N^Two), py.tuple({k, 1, N+1}));
soln = s.doit();
printing.print_latex(soln)

which give an answer of

$$ \frac{- \frac{\left(-1\right)^{N + 2}}{2} - \frac{1}{2}}{N^{2}} $$

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

$$ f(t) = e^{-a t}\cos{\omega t} $$

Python: 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

#!/usr/bin/env python
from sympy import exp, cos, laplace_transform, symbols
from sympy.printing.latex import print_latex
s = symbols('s')
t = symbols('t')
a = symbols('a', real=True, positive=True)
omega = symbols('omega', real=True)
f = exp(-a*t)*cos(omega*t)
L = laplace_transform(f, t, s)
print_latex(L[0])

MATLAB: 

1
2
3
4
5
6
7
8
9

sympy = py.importlib.import_module('sympy');
printing = py.importlib.import_module('sympy.printing.latex');
s = sympy.symbols('s');
t = sympy.symbols('t');
a = sympy.symbols('a', pyargs('real','True', 'positive','True'));
omega = sympy.symbols('omega', pyargs('real','True'));
f = sympy.exp(-a*t)*sympy.cos(omega*t);
L = sympy.laplace_transform(f, t, s); % a tuple, solution in 1st position
printing.print_latex(L{1})

which produces LaTeX that renders as

$$ \frac{a + s}{\omega^{2} + \left(a + s\right)^{2}} $$

Table of Contents | Previous: Solve simultaneous nonlinear equations Next : Accelerate MATLAB with Python and Numba

by Al Danial

Solve Systems of Nonlinear Equations in MATLAB with scipy.optimize.root()

Python MATLAB system of nonlinear equations GPS multivariate roots

Part 9 of the Python Is The Ultimate MATLAB Toolbox series.

Introduction

Systems of nonlinear equations appear in many diciplines of science, engineering, and economics1. Examples include mechanical linkages, power flow, heat transfer, and Calvo pricing. MATLAB’s Optimization Toolbox has functions that can solve such equations. If you don’t have access to this toolbox and you’re willing to include Python in your MATLAB workflow, you can use SciPy’s scipy.optimize.root() to find solutions to systems of nonlinear equations in MATLAB.

Example: The GPS Equations

Satellites in the global positioning system (GPS) constellation continuously transmit a navigation message which includes information about their orbital positions and the precise time. A GPS receiver that can get signals from at least four GPS satellites can compute its location by solving the nonlinear algebraic equations defining the intersection of four spheres, where the spheres' radii are determined by how long it took the messages to reach the receiver while traveling at the speed of light. Technically, only three spheres are needed to resolve two points in space, one of which can typically be rejected as a candidate based on its distance from the earth’s surface. However, since GPS receivers' clocks are inaccurate compared to the atomic clocks on the satellites, using a fourth satellite provides an extra equation that lets us additionally solve for the receiver’s clock error.

The four spheres problem can actually be simplified to a system of linear equations which is much easier to solve. For the sake of this post, however, the nonlinear equations are solved directly to demonstrate how such problems are set up in Python and MATLAB.

The equations for the four spheres are

$$ \begin{align*} (x-A_1)^2 + (y-B_1)^2 + (z-C_1)^2 &= (c(t_1-T_R+e))^2 \\ (x-A_2)^2 + (y-B_2)^2 + (z-C_2)^2 &= (c(t_2-T_R+e))^2 \\ (x-A_3)^2 + (y-B_3)^2 + (z-C_3)^2 &= (c(t_3-T_R+e))^2 \\ (x-A_4)^2 + (y-B_4)^2 + (z-C_4)^2 &= (c(t_4-T_R+e))^2 \end{align*} $$

where numeric subscripts denote the satellite; $A$, $B$, and $C$ are earth-centered coordinates (such as ECI or ECEF) of each satellite; $t$ is the time when the satellite sent its navigation message; $T_R$ is the receiver’s time; $e$ is the error in the receiver’s clock; and $c$ is the speed of light, 299,792.458 km/s. These values come from the GMU article:

satellite $A_i$ [km] $B_i$ [km] $C_i$ [km] $\Delta_i = t_i - T_R$ [seconds]
1 15600 7540 20140 0.07074
2 18760 2750 18610 0.07220
3 17610 14630 13480 0.07690
4 19170 610 18390 0.07242

We’re solving for the receiver’s position, ($x$, $y$, $z$), and the error in its clock, $e$. Our location ($x$, $y$, $z$) will be in the same coordinate system, ECEF or ECI, used to define ($A$, $B$, $C$) and $e$ will be in seconds.

Defining the function to solve

scipy.optimize.root()’s primary arguments are the function to solve and a starting guess. In this context, “solve” means to find values for the function’s inputs so that the function returns zeros for all outputs, in other words, to find the roots of the function. To this end we’ll rewrite our equations so they equal zeros then write a Python function that returns the four computed values of the left hand sides:

$$ \begin{align*} (x-A_1)^2 + (y-B_1)^2 + (z-C_1)^2 - (c(t_1-T_R+e))^2 &= 0\\ (x-A_2)^2 + (y-B_2)^2 + (z-C_2)^2 - (c(t_2-T_R+e))^2 &= 0\\ (x-A_3)^2 + (y-B_3)^2 + (z-C_3)^2 - (c(t_3-T_R+e))^2 &= 0\\ (x-A_4)^2 + (y-B_4)^2 + (z-C_4)^2 - (c(t_4-T_R+e))^2 &= 0 \end{align*} $$

The input variables $(x, y, z, e)$ need to be stored in a list or NumPy array. We’ll use the Python variable x as a four item array so that $(x, y, z, e)$ are stored in x[0] through x[3].

The function also needs each satellite’s position, $(A_i, B_i, C_i)$, and the time difference between the satellite’s clock and the receiver’s clock, $t_i - T_R$ which we’ll store as Dt:

Python 

# file:  gps_location.py
def F(x, A, B, C, Dt):
    c_light = 299792.458 # km/s
    return (
       (x[0]-A[0])**2 + (x[1]-B[0])**2 + (x[2]-C[0])**2 - (c_light*(Dt[0] + x[3]))**2,
       (x[0]-A[1])**2 + (x[1]-B[1])**2 + (x[2]-C[1])**2 - (c_light*(Dt[1] + x[3]))**2,
       (x[0]-A[2])**2 + (x[1]-B[2])**2 + (x[2]-C[2])**2 - (c_light*(Dt[2] + x[3]))**2,
       (x[0]-A[3])**2 + (x[1]-B[3])**2 + (x[2]-C[3])**2 - (c_light*(Dt[3] + x[3]))**2,
    )

We’ll also call this function from MATLAB so it needs to be in a file we can import. The examples below, both Python and MATLAB, will use F() from the file gps_location.py.

GPS Equations Solved in Python

We can solve for our position and clock error by populating arrays with satellite data then calling SciPy’s multivariate root finder with a starting guess. The earth’s radius is about 6378 km so a reasonable guess would be a random point near the earth’s surface, for example (6000, 0, 0). The clock error should be much less than a second so we’ll use 0 as a guess for that term.

The Python function we’re finding roots of, F(), takes four additional arguments, namely the three arrays of satellite position, and time deltas, in addition to the input vector x. These additional arguments can be passed through to F() by giving them to scipy.optimize.root() via the args keyword option:

Python 

#!/usr/bin/env python
import numpy as np
import scipy.optimize
from gps_function import F

A  = np.array([15600, 18760, 17610, 19170.0]) # km/s
B  = np.array([ 7540,  2750, 14630,   610.0]) # km/s
C  = np.array([20140, 18610, 13480, 18390.0]) # km/s
Dt = np.array([0.07074, 0.07220, 0.07690, 0.07242]) # s

guess = [6000, 0, 0, 0]

location = scipy.optimize.root(F, guess, args=(A, B, C, Dt))
print(location)

The output includes a text string indicating whether or not the solution converged (.message), the number of times the function was evaluated (.nfev), the $4 \times 4$ Jacobian matrix (.fjac) and the 4 terms of function at the computed solution (.fun). The solution itself is in the return variable’s .x attribute:

    fjac: array([[-0.42964552, -0.5286645 , -0.49262911, -0.54151189],
       [ 0.20220688, -0.32707893,  0.75636015, -0.52919852],
       [ 0.8684367 , -0.05711611, -0.42770859, -0.24417374],
       [-0.14260014,  0.78119842, -0.0479713 , -0.60588199]])
     fun: array([-1.19209290e-07, -5.96046448e-08, -1.19209290e-07, -5.96046448e-08])
 message: 'The solution converged.'
    nfev: 14
     qtf: array([ 3.35605879e-05, -1.62470220e-06, -2.25947450e-06,  2.73577863e-07])
       r: array([ 6.38261526e+04,  2.49011922e+04,  3.81328894e+04,  2.78876134e+10,
       -2.27583853e+04,  4.74674593e+03, -1.98003469e+09, -1.01001959e+04,
       -1.33002639e+09,  2.38095380e+08])
  status: 1
 success: True
       x: array([-4.17727096e+01, -1.67891941e+01,  6.37005956e+03,  3.20156583e-03])
(-1.1920928955078125e-07, -5.960464477539063e-08, -1.1920928955078125e-07, -5.960464477539063e-08)

From location.x we see that our GPS receiver is near the North Pole at (-41.772, -16.789, 6370.060) km and its clock is off by 3.20156583e-03 seconds. location.fun has the function’s value at the solution. Since all values should be zero, this array represents the solution’s error.

GPS Equations Solved in MATLAB

MATLAB can call scipy.optimize.root() as easily as Python can, provided that the function we’re finding the roots of is implemented in Python. (This is the same restriction we encountered minimizing a function in MATLAB with scipy.optimize.differential_evolution().) In our case we’ll reuse the function F() from gps_function.py:

MATLAB 2020b+ 

% file:  solve_gps.m
so  = py.importlib.import_module('scipy.optimize');
gps = py.importlib.import_module('gps_function');

A  = [15600, 18760, 17610, 19170.0]; % km/s
B  = [ 7540,  2750, 14630,   610.0]; % km/s
C  = [20140, 18610, 13480, 18390.0]; % km/s
Dt = [0.07074, 0.07220, 0.07690, 0.07242]; % s

guess = [6000, 0, 0, 0];
ABCDt = py.tuple({A, B, C, Dt});

location = so.root(gps.F, guess, pyargs('args',ABCDt))

The location found in MATLAB is the same as the one found in Python:

>> solve_gps

location =

  Python OptimizeResult with no properties.

        fjac: array([[-0.42964552, -0.5286645 , -0.49262911, -0.54151189],
           [ 0.20220688, -0.32707893,  0.75636015, -0.52919852],
           [ 0.8684367 , -0.05711611, -0.42770859, -0.24417374],
           [-0.14260014,  0.78119842, -0.0479713 , -0.60588199]])
         fun: array([-1.19209290e-07, -5.96046448e-08, -1.19209290e-07, -5.96046448e-08])
     message: 'The solution converged.'
        nfev: 14
         qtf: array([ 3.35605879e-05, -1.62470220e-06, -2.25947450e-06,  2.73577863e-07])
           r: array([ 6.38261526e+04,  2.49011922e+04,  3.81328894e+04,  2.78876134e+10,
           -2.27583853e+04,  4.74674593e+03, -1.98003469e+09, -1.01001959e+04,
           -1.33002639e+09,  2.38095380e+08])
      status: 1
     success: True
           x: array([-4.17727096e+01, -1.67891941e+01,  6.37005956e+03,  3.20156583e-03])

Note the object’s class: Python OptimizeResult. Extracting values from such an object is possible but clumsy as the return values are also Python objects:

>> format long e
>> location.get('x')

ans =

  Python ndarray:

    -4.177270957085472e+01    -1.678919410652270e+01     6.370059559223343e+03     3.201565829594195e-03

    Use details function to view the properties of the Python object.

    Use double function to convert to a MATLAB array.

We’d much prefer a native MATLAB struct. That’s easily done with py2mat.m 2:

>> loc = py2mat(location)

  struct with fields:

          x: [-4.177270957085472e+01 -1.678919410652270e+01 6.370059559223343e+03 3.201565829594195e-03]
    success: 1
     status: 1
       nfev: 14
       fjac: [4×4 double]
          r: [6.382615261849605e+04 2.490119219503466e+04 3.813288936837602e+04 2.788761343528719e+10 … ]
        qtf: [3.356058790625163e-05 -1.624702195475986e-06 -2.259474496258976e-06 2.735778627939572e-07]
        fun: [-1.192092895507812e-07 -5.960464477539062e-08 -1.192092895507812e-07 -5.960464477539062e-08]
    message: "The solution converged."

Table of Contents | Previous: Solve global optimization problems | Next : Symbolic math

  1. Wolfgang Eichhorn & Winfried Gleißner, 2016. “Nonlinear Functions of Interest to Economics. Systems of Nonlinear Equations,” Springer Texts in Business and Economics, in: Mathematics and Methodology for Economics, edition 1, chapter 7, pages 301-371, Springer. ↩︎

  2. I added support for Python OptimizeResult objects in py2mat.m on 2022-09-16. Earlier versions will fail to convert such objects. ↩︎