SymPy

Scientific Computing

Author

Prof. Calvin

Why SciPy?

What is SymPy?

A Python library for symbolic mathematics. It aims to become a full-featured computer algebra system (CAS) while keeping the code as simple as possible in order to be comprehensible and easily extensible. SymPy is written entirely in Python.

  • Symbolic Python

Why SymPy

  • “SymPy is free both as in speech and as in beer.”
  • “Most computer algebra systems invent their own language. Not SymPy.”
  • “An advantage of SymPy is that it is lightweight.”

Why not SymPy

  • To my knowledge, there’s no real competitors with SymPy.
  • The closest are probably SageMath and Mathematica.
  • I use browser-based Mathematica via https://www.wolframalpha.com/ from time to time (when I don’t have a Python installation handy).

Cite

@article{10.7717/peerj-cs.103,
     title = {SymPy: symbolic computing in Python},
     author = {Meurer, Aaron and Smith, Christopher P. and Paprocki, Mateusz and \v{C}ert\'{i}k, Ond\v{r}ej and Kirpichev, Sergey B. and Rocklin, Matthew and Kumar, AMiT and Ivanov, Sergiu and Moore, Jason K. and Singh, Sartaj and Rathnayake, Thilina and Vig, Sean and Granger, Brian E. and Muller, Richard P. and Bonazzi, Francesco and Gupta, Harsh and Vats, Shivam and Johansson, Fredrik and Pedregosa, Fabian and Curry, Matthew J. and Terrel, Andy R. and Rou\v{c}ka, \v{S}t\v{e}p\'{a}n and Saboo, Ashutosh and Fernando, Isuru and Kulal, Sumith and Cimrman, Robert and Scopatz, Anthony},
     year = 2017,
     month = jan,
     keywords = {Python, Computer algebra system, Symbolics},
     abstract = {
                SymPy is an open source computer algebra system written in pure Python. It is built with a focus on extensibility and ease of use, through both interactive and programmatic applications. These characteristics have led SymPy to become a popular symbolic library for the scientific Python ecosystem. This paper presents the architecture of SymPy, a description of its features, and a discussion of select submodules. The supplementary material provide additional examples and further outline details of the architecture and features of SymPy.
             },
     volume = 3,
     pages = {e103},
     journal = {PeerJ Computer Science},
     issn = {2376-5992},
     url = {https://doi.org/10.7717/peerj-cs.103},
     doi = {10.7717/peerj-cs.103}
    }

pip again

  • Just like NumPy, Matplotlib is a Python package which we install via pip
python3 -m pip install sympy
  • That might take a moment, when it does we can check it worked!

Motivation

Example

  • Computers don’t hold numbers that precisely.
10.0 ** 100 == 10.0 ** 100 + 47
True

Fractions

  • Much worse with fractions I find.
1 / 7000000 * 7000000

\(\displaystyle 1.0\)

Floats:

  • The IEEE 754 Floating-Point Standard

  • It is basically scientific notation that fits in a fixed amount of characters.

What is it?

  • The IEEE 754 standard defines formats for representing floating-point numbers.
  • It specifies how floating-point numbers are stored and operated on in computer hardware.
  • Most modern CPUs adhere to this standard.
import sys
# Check the floating-point precision on your system
sys.float_info
sys.float_info(max=1.7976931348623157e+308, max_exp=1024, max_10_exp=308, min=2.2250738585072014e-308, min_exp=-1021, min_10_exp=-307, dig=15, mant_dig=53, epsilon=2.220446049250313e-16, radix=2, rounds=1)

Why is it important?

  • Ensures portability and consistency of numerical computations across different systems.
  • Without a standard, the same calculation could yield different results on different machines.
  • Essential for reliable scientific and engineering software.

Single Precision (Float32)

  • Uses 32 bits to represent a number.
  • 1 sign bit, 8 exponent bits, 23 significand (mantissa) bits.
  • Offers approximately 7 decimal digits of precision.
import numpy as np
np.float32(1/7)
np.float32(0.14285715)

Double Precision (Float64)

  • Uses 64 bits to represent a number.
  • 1 sign bit, 11 exponent bits, 52 significand bits.
  • Offers approximately 15-17 decimal digits of precision. This is the default in Python and NumPy.
# Example of a float64 number (default in Python)
1/7

\(\displaystyle 0.142857142857143\)

Special Values

  • Infinity \(\infty\): Result of overflow or division by zero.
lil = np.finfo(np.float64).resolution # smallest recognizable value
big = np.finfo(np.float64).max # biggest recognizable value
big / lil
C:\Users\cd-desk\AppData\Local\Temp\ipykernel_21564\563832135.py:3: RuntimeWarning:

overflow encountered in scalar divide

\(\displaystyle \infty\)

  • Not a Number (NaN): Result of undefined operations (e.g., \(0/0\), \(\sqrt{-1}\)).
np.sqrt(-1)
C:\Users\cd-desk\AppData\Local\Temp\ipykernel_21564\3438155168.py:1: RuntimeWarning:

invalid value encountered in sqrt

\(\displaystyle \text{NaN}\)

Pitfall 1: Limited Precision

  • Not all real numbers can be represented exactly.
  • Decimal numbers like 0.1 often have an infinitely repeating binary representation.
  • This leads to small, unavoidable rounding errors.
.1 + .2

\(\displaystyle 0.3\)

Pitfall 2: Accumulation

  • Small rounding errors can accumulate over many operations.
  • This can lead to significant inaccuracies in long computations or iterative algorithms.
  • Careful algorithm design and error analysis are crucial.
total = 0.0
for _ in range(100000):
    total += 0.1
total, 100000 * 0.1

\(\displaystyle \left( 10000.0000000188, \ 10000.0\right)\)

Pitfall 3: Comparison Issues

  • Due to limited precision, direct equality comparisons (==) between floating-point numbers are often unreliable.
(1 + 2) / 10

\(\displaystyle 0.3\)

1 / 10 + 2 / 10

\(\displaystyle 0.3\)

Pitfall 4: Subtractive Cancellation

  • Occurs when subtracting two nearly equal numbers.
  • The most significant bits cancel out, leaving only the less significant bits.
  • Can drastically reduce the effective precision.
x = 1.0000000000000001
y = 1.0000000000000000
x - y

\(\displaystyle 0.0\)

Symbolic Computation

Credit

Quoth SymPy

Symbolic computation deals with the computation of mathematical objects symbolically. This means that the mathematical objects are represented exactly, not approximately, and mathematical expressions with unevaluated variables are left in symbolic form.

Non-symbolic Computation

  • A common problem in computing is to compute the likes of distance between points separated by vectors.
  • Imagine two objects are separated by a horizontal displacement of 3 units, and vertical displacement of 4 units, and we wish to determine the minimum distance.
  • A straightforward application of the Pythagorean Theorem.
np.sqrt(3 * 3 + 4 * 4)

\(\displaystyle 5.0\)

Pitfalls

  • These distances quickly become inaccurate.
np.sqrt(2 * 2 + 4 * 4)

\(\displaystyle 4.47213595499958\)

  • Then stuff like this happens:

Symbolic Computation

  • In point of fact, the solution to \(\sqrt{2^2 + 4*2}\) is actually \(\sqrt{20}\) and there’s no other graceful way to represent it.
  • So we’ll use symbols to store values.
import sympy
x = sympy.symbols('x')
y = sympy.symbols('y')
dist = sympy.sqrt(x + y)
dist

\(\displaystyle \sqrt{x + y}\)

Expressions

  • We term this a “symbolic expression”.
  • We can perform operations…
2 * dist

\(\displaystyle 2 \sqrt{x + y}\)

  • We can add other symbols…
z = sympy.symbols('z')
z + dist

\(\displaystyle z + \sqrt{x + y}\)

Niceties

  • We can declare multiple symbols at once:
a, b, c = sympy.symbols('a b c')
  • SymPy will automatically simplify.
sympy.sqrt(20)

\(\displaystyle 2 \sqrt{5}\)

  • Compare:
np.sqrt(20), np.sqrt(5), 2 * np.sqrt(5)

\(\displaystyle \left( 4.47213595499958, \ 2.23606797749979, \ 4.47213595499958\right)\)

Operations

Extended Floats

  • SymPy uses floats by default, but doesn’t have to.
  • You can specify a number of decimals for any value.
  • We also should that SymPy contains some useful constants!
# "Evaluate as float to 100 digits"
# We use sympy.py, not e.g. np.pi
sympy.pi.evalf(100)

\(\displaystyle 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068\)

Substitution

  • Often we want to solve an equation algebraically and also know a numerical solution.
  • The usefulness of SymPy is to do both, and maximally simply the result to minimize error.
  • We can use subs() to get solutions given values.
dist

\(\displaystyle \sqrt{x + y}\)

Value Substitution

  • Substitute 3 for x
dist.subs(x, 3)

\(\displaystyle \sqrt{y + 3}\)

Expression Substitution

  • Perhaps y is, itself, a distance expressed over a right triangle with sides a and b
dist.subs(y,a*a+b*b)

\(\displaystyle \sqrt{a^{2} + b^{2} + x}\)

Strings

  • Sometimes we want to take a Python expression and convert to a symbolic SymPy expression.
str_expr = "x**2 + 3*x - 1/2"
from sympy import sympify
expr = sympify(str_expr)
expr

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

  • Then calculate with .subs()
expr.subs(x, 2)

\(\displaystyle \frac{19}{2}\)

Functions

  • Sibling of vectorize
  • Takes a SymPy expression, makes a Python function.
expr # so we remember.

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

  • Calculate…
f = sympy.lambdify(x, expr, "scipy")
f(np.arange(10))
array([ -0.5,   3.5,   9.5,  17.5,  27.5,  39.5,  53.5,  69.5,  87.5,
       107.5])

Aside: SciPy/NumPy

  • Always use lambdify specified with "scipy" (or "numpy") if you have SciPy or even just NumPy installed
  • It uses NumPy’s more powerful (than Python’s) mathematical operations.

Aside: Plots

# import sympy
# from sympy.abc import x
# import scipy
import matplotlib.pyplot as plt
f = sympy.lambdify(x,x*x,"scipy")
xs = np.linspace(-5,5) 
plt.plot(xs, f(xs))

Simplification

  • If I use SymPy, I don’t have to remember this:
sympy.simplify(sympy.sin(x)**2 + sympy.cos(x)**2)

\(\displaystyle 1\)

  • Or figure out this:
sympy.simplify((x**3 + x**2 - x - 1)/(x**2 + 2*x + 1))

\(\displaystyle x - 1\)

Polynomials

  • Often we don’t want a simplified polynomial.
sympy.simplify(x**2 + 2*x + 1)

\(\displaystyle x^{2} + 2 x + 1\)

  • We perhaps instead wish to factor()
sympy.factor(x**2 + 2*x + 1)

\(\displaystyle \left(x + 1\right)^{2}\)

Expand

  • We may also start with factored form.
eq = x**2 - (x+1) * (x-1)
  • This can deceptively by made smaller by expanding it.
sympy.expand(eq)

\(\displaystyle 1\)

Fractions

  • We can also cancel out fractions.
eq = (x**2 + 2*x + 1)/(x**2 + x)
eq

\(\displaystyle \frac{x^{2} + 2 x + 1}{x^{2} + x}\)

  • Using sympy.factor()
sympy.factor(eq)

\(\displaystyle \frac{x + 1}{x}\)

Powers

  • SymPy can also perform some simplifications related to exponention.
  • Movie size in memory is frame size time frame rate times duration.
eq = x**10*x**9*x**5*x*9 # Avatar in HD is ~4 hrs long at ~60fps
eq

\(\displaystyle 9 x^{25}\)

  • Or just sub in 2
eq.subs(x,2)

\(\displaystyle 301989888\)

Calculus

Derivatives

  • We can calculate derivatives.
sympy.diff(x*x, x)

\(\displaystyle 2 x\)

Integrals

sympy.integrate(sympy.sin(x), x)

\(\displaystyle - \cos{\left(x \right)}\)

Evaluate Integrals

  • Give a variable a domain…
sympy.integrate(sympy.exp(x), (x, 0, 1))

\(\displaystyle -1 + e\)

Limits

  • One famous problem is to compute the value of the following function as it approaches zero: \[ \lim_{x\to 0} \frac{\sin{x}}{x} \]

Evaluate limits

eq = sympy.sin(x)/x
eq

\(\displaystyle \frac{\sin{\left(x \right)}}{x}\)

  • We provide a variable and the value it approaches.
sympy.limit(eq, x, 0)

\(\displaystyle 1\)

Viewing limits

  • To create a printable limit object that can be evaluated.
    • Capitalize
    • Resolve via .doit
    • Also works for Derivative and Integral
lm = sympy.Limit(eq, x, 0)
lm, lm.doit()

\(\displaystyle \left( \lim_{x \to 0^+}\left(\frac{\sin{\left(x \right)}}{x}\right), \ 1\right)\)

Solutions

Solve Equations

The Python package SymPy can symbolically solve equations, differential equations, linear equations, nonlinear equations, matrix problems, inequalities, Diophantine equations, and evaluate integrals. SymPy can also solve numerically.

Algebra

  • We use the much nicer way to get variable names.
from sympy.abc import x, y
  • We will aim to solve the following: \[ x^2 = y \]

Write in SymPy

  • SymPy, like many algebra systems, expects equations to be expressed as equal to zero.
  • We note the following equivalence. \[ x^2 = y \equiv x^2 - y = 0 \]
  • So we write:
eq = x*x - y
eq

\(\displaystyle x^{2} - y\)

Solve

  • To solve, we simply use sympy.solve()
    • We specify the equation to solve, and
    • The variable for which to solve.
sympy.solve(eq, x)

\(\displaystyle \left[ - \sqrt{y}, \ \sqrt{y}\right]\)

  • We could also do:
sympy.solve(eq, y)

\(\displaystyle \left[ x^{2}\right]\)

== and Eq()

  • SymPy, for a variety of good reasons, struggles with == in equations, so use Eq()
eq = sympy.Eq(x*x,y)
eq

\(\displaystyle x^{2} = y\)

  • It works as expected.
sympy.solve(eq)

\(\displaystyle \left[ \left\{ y : x^{2}\right\}\right]\)

Aside: Rationals

  • The is true for fractions with Rational()
sympy.Rational(5/2)

\(\displaystyle \frac{5}{2}\)

  • In general, if something would be “weird” with Sympy, there is a capitalized name that does what you would expect.
  • Just check the documentation.

Printing

Pretty Print

  • SymPy provides a number of ways to see equations.
  • By default, these slides use “LaTeX”, which I regard as the standard for mathematical typesetting.
  • However, it doesn’t often show up that way in the terminal!

init_printing

  • Here’s what my terminal looks like:
>>> import sympy
>>> from sympy.abc import x, y, z
>>> sympy.Eq(x*x,y)
Eq(x**2, y)
>>> sympy.init_printing()
>>> sympy.Eq(x*x,y)
 2
x  = y

Options

  • The following visualization methods are usably by SymPy:
    • str
    • srepr
    • ASCII pretty printer
    • Unicode pretty printer
    • LaTeX
    • MathML
    • Dot

Example

  • We’ll use an integral example:
sympy.Integral(sympy.sqrt(1/x), x)

\(\displaystyle \int \sqrt{\frac{1}{x}}\, dx\)

str

  • Just gives the string, basically as we typed it in.
str(sympy.Integral(sympy.sqrt(1/x), x))
'Integral(sqrt(1/x), x)'

srepr

  • srepr is string representation and is more exact and verbose.
  • I don’t use it often, but it can be helpful to examine equations when I get unexpected answers.
sympy.srepr(sympy.Integral(sympy.sqrt(1/x), x))
"Integral(Pow(Pow(Symbol('x'), Integer(-1)), Rational(1, 2)), Tuple(Symbol('x')))"

ASCII

  • ASCII is the oldest major standard for character displays, with 127 characters including non-printing or whitespace characters like tab and space.
  • We can make “ASCII art” of equations via pprint.
    • We must set use_unicode to False
sympy.pprint(sympy.Integral(sympy.sqrt(1/x), x), use_unicode=False)
  /          
 |           
 |     ___   
 |    / 1    
 |   /  -  dx
 | \/   x    
 |           
/

Unicode

  • We get some slightly smoother lines with the more complete Unicode character set.
    • Unicode is the modern more complete character set including non-English characters and the likes of emojis. 🤔💭🔢✖️🧮
sympy.pprint(sympy.Integral(sympy.sqrt(1/x), x), use_unicode=True)
⌠           
⎮     ___   
⎮    ╱ 1    
⎮   ╱  ─  dx
⎮ ╲╱   x    
⌡

LaTeX

  • My favorite (by far!)
  • I use LaTeX (.tex) a lot in these and other slides.
sympy.latex(sympy.Integral(sympy.sqrt(1/x), x))
'\\int \\sqrt{\\frac{1}{x}}\\, dx'
  • I don’t see the double backslash in my terminal, that is just an artiffact of the slides.
>>> sympy.latex(sympy.Integral(sympy.sqrt(1/x), x))
\int \sqrt{\frac{1}{x}}\, dxyy

Using LaTeX

  • I can directly display that LaTeX in these slides…
  • I type:
file.tex
$$
\int \sqrt{\frac{1}{x}}\, dxyy
$$
  • $$ denotes LaTeX “math mode”, versus just basic text editing.
  • We see:

\[ \int \sqrt{\frac{1}{x}}\, dxyy \]

QuickLaTeX

\int \sqrt{\frac{1}{x}}\, dxyy
  • I see this:

Overleaf

tex command

  • I write the LaTeX output to file.
integral.py
import sympy 
x = sympy.symbols('x')
print(sympy.latex(sympy.Integral(sympy.sqrt(1/x), x)))
  • And then:
$ python3 integral.py > file.tex
  • Or all at once using ; to separate lines
$ python3 -c "import sympy; x = sympy.symbols('x'); print(sympy.latex(sympy.Integral(sympy.sqrt(1/x), x)))" > file.tex

nvim edits

  • I open the file with nvim and add lines specifying I want a document and that I want “math mode”
  • Before:
\int \sqrt{\frac{1}{x}}\, dx
  • After:
\begin{document}
$$
\int \sqrt{\frac{1}{x}}\, dx
$$
\end{document}

Render

  • Simply render to a .pdf via
$ pdflatex file.tex

ODEs

Ordinary Differential Equations

A differential equation is an equation that relates one or more unknown functions and their derivatives.

An ordinary differential equation (ODE) is a differential equation (DE) dependent on only a single independent variable.

Credit

Equation

  • Consider the decay of tritium as an example.

\[ {}^3\text{H} \xrightarrow{\lambda} {}^3\text{He} + \text{e}^- + \bar{\nu}_\text{e} \]

  • This is, by the way, LaTeX:
tritium.tex
{}^3\text{H} \xrightarrow{\lambda} {}^3\text{He} + \text{e}^- + \bar{\nu}_\text{e}

Tritium

Tritium is a radioactive isotope of hydrogen. It has the same number of protons and electrons as hydrogen but has 2 neutrons, whereas regular hydrogen does not have any. This makes tritium unstable and radioactive. Tritium is produced naturally from interactions of cosmic rays with gases in the upper atmosphere, and is also a by-product of nuclear reactors.

Tritium

Like all radioactive isotopes, tritium decays. As it decays, it emits beta radiation.

The physical half-life of tritium is 12.33 years, meaning that it takes just over 12 years for tritium to decay to half of its original amount. As tritium decays, it changes to helium.

Derivatives

  • We denote the “number density” (number of atoms in a sample) of \(^3\text{H}\) as a function of time \(y(t)\)
  • The rate of change (derivate)

\[ \frac{d}{d t} y{\left(t \right)} = - \lambda y{\left(t \right)} \]

  • This is a differential equation because an equation is a function of its own derivative.

Exercise

  1. Specify this equation in SymPy
  2. Solve useing sympy.dsolve()

Tools

  • This is how I got lambda:
t, l = sympy.symbols('t lambda') # use "l" as "lambda" is a Python keyword
  • This is \(y(t)\)
y = sympy.Function('y')(t) # use capital Function ala Rational, Integral
  • Take derivatives as follows:
f = sympy.Function('f')(x)
f.diff(x)

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

Solution

  • Write:
Code 
t, l = sympy.symbols('t lambda')
y = sympy.Function('y')(t)
dydt = y.diff(t)
expr = sympy.Eq(dydt, -l*y)
expr

\(\displaystyle \frac{d}{d t} y{\left(t \right)} = - \lambda y{\left(t \right)}\)

  • Solve:
Code 
sympy.dsolve(expr)

\(\displaystyle y{\left(t \right)} = C_{1} e^{- \lambda t}\)