by Al Danial

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

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

Join me again on October 29, 2022 to see how Python with Numba can make MATLAB programs run much faster.