Lesson A1: Symbolic Mathematics (1)

Instructor: Yuki Oyama, Prprnya
The Christian F. Weichman Department of Chemistry, Lastoria Royal College of Science

This material is licensed under CC BY-NC-SA 4.0

Welcome to the second part of the course! The first four lessons covered the basics of Python and its core libraries, so they’re meant to be learned in sequence. This second part, however, is more application-oriented, and the lessons are independent of each other. We’ve arranged them alphabetically—not as a learning order, but simply to keep track of the order in which they were created. From here on, you are free to choose any lesson to start with( ´▽｀)

Introduction¶

Perhaps all calculations we’ve done in previous lessons have been done using numerical values. However, in many cases, we want to do calculations using symbols. For example, given the Maxwell–Boltzmann distribution:

f(v) = 4\pi \left( \frac{m}{2\pi k_B T} \right)^\frac{3}{2} v^2 e^{-\frac{mv^2}{2k_BT}}

(1)

with $v$ the speed of a particle, $m$ the mass of this particle, $T$ the temperature, and $k_B$ the Boltzmann constant. We want to find the most probable speed $v_\text{mp}$ , which is the speed at which the particle is most likely to be at a given temperature. The canonical way to do this is to find the maximum of $f(v)$ , and this requires us to evaluate the derivative of $f(v)$ with respect to $v$ , and set it to be zero. What if we don’t want to do this tedious (well, it’s not quite tedious—you should know it! but I’m lazy...) calculation? You can use an expensive calculator with CAS (Computer Algebra System) functionality to do this, or you can use software like Mathematica or Maple... These all take a lot of money! Luckily, there is a fantastic library in Python: SymPy, that can carry symbolic mathematics like finding derivatives, and this is what we will study in this lesson.

By the way, the most probable speed is $v_\text{mp} = \sqrt{\frac{2k_BT}{m}}$ . Do you get it?

Exact is Magic!¶

What an amaze symbolic mathematics is it can give us exact answers, not just numbers, to problems we’d otherwise have to work out by hand. Suppose you want to get the simplified form of $\sqrt{8}$ (this time... you MUST know this!). Let’s see what NumPy tells us…

import numpy as np
np.sqrt(8)

It’s a number… well, kinda not the stuff we want. Let’s try SymPy... oh, we need to import it first:

import sympy as sp

sp is the recommended alias for sympy.

This time let’s see:

sp.sqrt(8)

Wow! This is an exact result! We can even do some more complicated things:

\sqrt{8} + \sqrt{2}

(2)

sp.sqrt(8) + sp.sqrt(2)

... even this:

\int \left[e^\xi \sin(\xi) + e^\xi \cos(\xi)\right] \, d\xi

(3)

xi = sp.symbols('xi')
sp.integrate(sp.exp(xi)*sp.sin(xi) + sp.exp(xi)*sp.cos(xi), xi)

Also, do you notice that the results are beautiful mathematical expressions, just like the LaTeX expressions we’ve seen before? This is because SymPy has this feature called “pretty printing,” which automatically renders the output to make it look nice.

Symbols and Functions¶

Symbols¶

Okay, fancy demonstrations over. Let’s get back to the real deal: symbols. Have you noticed that when we are going to carry out calculations of algebra—like $\xi$ above, we need to create a symbol using sp.symbols()? In SymPy, we cannot use x or y directly, but we need to assign them to be symbols. If you try to add a defined symbol (like xi above) with an undefined symbol, SymPy will throw an error:

xi + x

We need to define x as symbols first:

x = sp.symbols('x')
xi + x

You can assign multiple symbols at once:

x, y, z = sp.symbols('x y z')
x + y + z

The name of a symbol and the name of the variable it is assigned to need NOT have anything to do with one another. For example, we can assign a to be the symbol $b$ , and b to be the symbol $a$ :

a, b = sp.symbols('b a')
a # this should gives $b$

b # this should gives $a$

But! This is really a bad practice! Please don’t do this...

There is a shortcut to define common symbols from alphabet and Greek letters. That is, you can import sympy.abc to get the entire alphabetical (Roman) letters and lowercase Greek letters, except letters you have defined before (x, y, z, xi, a, b in this document) and some special letters reserved for Python (like lambda)

from sympy.abc import *

Here, import * means that all elements in sympy.abc will be imported. Now you can check:

(c + d) / e + (P + Q) + alpha * beta / gamma

Uppercase Greek letters need to be defined by our own:

Gamma = sp.symbols('Gamma')
Gamma

We can also assign specific predicates to symbols, like positive and negative. This can be achieved by giving assumptions when defining the symbol:

x1 = sp.symbols('x', positive=True)
x2 = sp.symbols('x', positive=True)
x1 == x2

See? x1 and x2 are not equal, because x1 is defined as a positive symbol, but x2 has no restrictions on it. Some commonly seen predicates are shown below.

Predicate	Definition	Implications / Relations
complex	A complex number $a + bi$ , where $a, b \in \mathbb{R}$ . All complex numbers are finite and include all real numbers.	→ `commutative`, → `finite`
real	A real number ( $\mathbb{R}$ ). Every real is also complex ( $\mathbb{R} \subset \mathbb{C}$ ); includes all rationals and integers.	→ `complex`, == (`negative`
imaginary	A number of the form $bi$ , where $b \in \mathbb{R}$ and $b \neq 0$ . Complex but not real.	→ `complex`, → `!real`
integer	An integer $(\cdots, -2, -1, 0, 1, 2, \cdots)$ .	→ `rational`, → `real`
even	An even integer ( $2n$ ). Includes zero.	→ `integer`, → `!odd`
odd	An odd integer ( $2n + 1$ ).	→ `integer`, → `!even`
prime	A positive integer greater than 1 that has no divisors other than 1 and itself.	→ (`integer` & `positive`)
nonzero	A real or complex number that is not zero.	== (`!zero`), → (`real`
positive	A real number $> 0$ . All positive numbers are finite.	== (`nonnegative` & `nonzero`), → `real`
negative	A real number $< 0$ . All negative numbers are finite.	== (`nonpositive` & `nonzero`), → `real`

A complete list of predicates can be accessed in the SymPy documentation of assumptions.

Functions¶

Functions can be defined by using Function class, and the syntax is similar to defining symbols. For example, we can define an abstract function $f(x)$ as follows:

f = sp.Function('f')
f

Well... this is only a function definition, not a function call. For calling a function, we need to use the f(x):

f(x)

Explicitly including variables can give the same result as above:

g = sp.Function('g')(x)
g

You can also define a function with multiple arguments:

h = sp.Function('h')(x, y)
h

Assumptions can be provided to a Function in the same way as they are for a Symbol. Alternatively, you can define a Symbol that already has assumptions, and when used as the function’s name, the function automatically inherits those assumptions.

f_real = sp.Function('f', real=True)
f_real(x).is_real

f_real_inherit = sp.Function(sp.Symbol('f', real=True))
f_real_inherit(x).is_real

You can also define a specific function in SymPy, but this is more complicated than defining a general function. We will not spend much time on this, but you can refer to the SymPy documentation of Writing Custom Functions.

Basic Operations¶

Substitution and Evaluation¶

One of the most important thing that you may want to do with symbols is to substitute them with numbers or other symbols. For example, recall that the wavefunction of a quantum harmonic oscillator is given by

\psi_n(x) = \left(\frac{m \omega}{\pi \hbar}\right)^{1/4} \frac{1}{\sqrt{2^n n!}} H_n\left(\sqrt{\frac{m \omega}{\hbar}} x\right) e^{-\frac{m \omega x^2}{2 \hbar}}

(4)

Sometimes people may shorten this writing by using $\xi = \sqrt{\frac{m \omega}{\hbar}} x$ , so the above expression reduces to

\psi_n(\xi) = \left(\frac{m \omega}{\pi \hbar}\right)^{1/4} \frac{1}{\sqrt{2^n n!}} H_n\left(\xi\right) e^{\frac{-\xi^2}{2}}

(5)

However, if we want to do the above procedure in reverse, I bet you wouldn’t want to manually substitute all the symbols twice and then simplify the square root of a square by hand. Using SymPy, we can do this easily. First, let’s define the expression:

# Don't be worried about these imports!
# They are just for defining constants and functions from sympy
# We will show you how to use them later
from sympy.functions.special.polynomials import hermite
from sympy.physics.quantum import hbar

psi_n = ((m * omega)/(sp.pi * hbar))**(1/4) * 1/sp.sqrt(2**n * sp.factorial(n)) * hermite(n, xi) * sp.exp(-xi**2 / 2)
psi_n

We can substitute the symbol $\xi$ with $\sqrt{\frac{m \omega}{\hbar}} x$ with the subs() function, which is a method, so you need to call it by add .subs() after your expression. This function takes two arguments: the symbol to be substituted, and the value to be substituted with. You can check the documentation of subs() for more details.

psi_n.subs(xi, sp.sqrt(m*omega/hbar) * x)

The result is a symbolic expression, which can be evaluated using evalf():

You can try substituiton again with $n=0$ to find the ground state wavefunction:

psi_0 = psi_n_full.subs(n, 0)
psi_0

It’s not difficult to generate excited state wavefunctions with different n:

from ipywidgets import interact, BoundedIntText

# Generate a widget to let user input n from 0 to 15
@interact(nval=BoundedIntText(value=0, description='n:', min=0, max=15))
def psi_n_full_expr(nval):
    return psi_n_full.subs(n, nval)

What if we want to evaluate the wavefunction $\psi_0$ at a specific point, like $\psi_0(0)$ ? We can use subs() again:

psi_n_full.subs([(x, 0), (n, 0)])

This will give us an exact result. We can also convert this expression to a numerical value, by using evalf():

# Normally we don't use `print()` to output SymPy results
# Here is a special case - the output needs to be converted to string by `print()`
# You can try not to use `print()` and see what happens
print(psi_n_full.subs([(x, 0), (n, 0)]).evalf())

Hmm... it seems that we forget to substitute $\omega$ and $m$ with their actual values. Assume that we are treating a $\ce{H2}$ molecule as an harmonic oscillator, we can substitute them with $\omega = 8.28 \times 10^{14}\, \mathrm{s^{-1}}$ and $m = 8.37 \times 10^{28}\, \mathrm{kg}$ :

psi_n_full.subs([(x, 0),
                 (n, 0),
                 (omega, 8.28e-14),
                 (m, 8.37e-28)]).evalf()

SymPy can evaluate floating point expressions to arbitrary precision. By default, 15 digits of precision are used, but you can pass any number as the argument to evalf(). For example:

psi_n_full.subs([(x, 0),
                 (n, 0),
                 (omega, 8.28e-14),
                 (m, 8.37e-28)]).evalf(5)

See how it gives the answer? Basically, you can use SymPy as a really powerful calculator( ´▽｀)

For more information about subs() and evalf(), you can refer to the SymPy tutorials of Basic Operations.

“Lambdification”¶

subs() and evalf() work well for quick evaluations, but they become inefficient when you need to compute values at many points. If you’re evaluating an expression thousands of times—especially at machine precision—SymPy will be much slower than necessary. In such cases, it’s better to use NumPy or SciPy instead.

The easiest way to make a SymPy expression numerically evaluable is with sp.lambdify(), which converts symbolic expressions into fast, NumPy-compatible functions. This function takes two arguments: the variables to be substituted, and the expression to be converted. If you want to convert the expression to a function that takes multiple arguments, just pass a tuple of variables to the first argument of sp.lambdify(). Let’s try it out with our harmonic oscillator wavefunction:

# We first define a `psi_for_lambdify` that substitutes the values of `omega` and `m`
psi_for_lambdify = psi_n_full.subs([(omega, 8.28e-14), (m, 8.37e-28)])
# Then we use `lambdify()` to convert it to a function that takes `x` and `n` as arguments
harmonic_psi_from_sympy = sp.lambdify((x, n), psi_for_lambdify)
# Finally, we can evaluate the function at any `x` and `n` we want
harmonic_psi_from_sympy(1, 0)

From here, we can pass an NumPy array to this function and we can easily plot the wavefunction with Matplotlib:

import matplotlib.pyplot as plt

x_vals = np.linspace(-10, 10, 100)
psi_vals = harmonic_psi_from_sympy(x_vals, 0)

plt.plot(x_vals, psi_vals)
plt.xlabel(r'$x$')
plt.ylabel(r'$\psi(x)$')
plt.show()

Nice plots, right? This is similar to what we did in the End-of-Lesson Problems of Lesson 2.

There are more functionalities of sp.lambdify(), like specifying the module to use for the numeric library, or the module to use for the printing library. You can check more details in the SymPy documentation of Lambdify.

Exercise : Use SymPy to make a function of our old friend—Morse potential:

V_\text{Morse}(r) = D_e \left( 1 - e^{-a \left( r - r_0 \right)} \right)^2

(6)

Evaluate it at $r = r_0$ , and use sp.lambdify() to make a function that takes $r$ as an argument and returns the value of $V_\text{Morse}(r)$ . Try to plot ths function with $D_e = 10.98\,\mathrm{eV}$ , $a = 2.32\,\mathrm{Å}^{-1}$ , and $r_0 = 1.128\,\mathrm{Å}$ from $r = 0.8\, \mathrm{Å}$ to $r = 3\, \mathrm{Å}$ .

Simplification¶

While “simplification” sounds like a great word, it actually covers a lot. We all know that expressions like $1 + 2x - x$ can be simplified to $1 + x$ , and $\cos^2 x + \sin^2 x = 1$ , but what about expressions like $x^2 + 2x + 1$ ? Do you think that’s already simplified? Well, if you want to get $(x + 1)^2$ , that’s factorization, but if you prefer to keep it as it is, that’s also a simplified—or rather, expanded form.

Naive Simplification¶

SymPy has many built-in simplification functions, each useful for different levels or types of simplification. The most common one is probably sp.simplify(). As its name suggests, this function simplifies expressions—and it does so smartly, using a heuristic approach. We’ll go over the details later, but first, let’s look at a few examples:

sp.simplify(sp.sin(x)**2 + sp.cos(x)**2)

\text{simplify}\; \sin^2 x + \cos^2 x

(7)

sp.simplify((x**3 + x**2 - x - 1)/(x**2 + 2*x + 1))

\text{simplify}\; \frac{x^3 + x^2 - x - 1}{x^2 + 2x + 1}

(8)

sp.simplify(sp.gamma(x)/sp.gamma(x - 2))

\text{simplify}\; \frac{\Gamma(x)}{\Gamma(x - 2)}

(9)

(Here $\Gamma(x)$ is the gamma function. It is closely related to factorials.)

What about $x^2 + 2x + 1$ ? In fact, sp.simplify() doesn’t change it—it either doesn’t know a simpler form or considers it already simplified, so it just returns the original expression:

sp.simplify(x**2 + 2*x + 1)

\text{simplify}\; x^2 + 2x + 1

(10)

Expansion¶

For polynomials like $x^2 + 2x + 1$ , SymPy has functions to simplify them specifically. sp.expand() will convert a polynomial to a canonical form of a sum of monomials. Well, $x^2 + 2x + 1$ is already in this form, so let’s see what happens if we use sp.expand() on $(x+1)^2$ :

sp.expand((x+1)**2)

\text{expand}\; (x+1)^2

(11)

sp.expand() might not sound like a simplification function. After all, its name suggests making expressions larger, not smaller, and that’s usually true, but sometimes applying sp.expand() can actually make an expression simpler due to term cancellations.

sp.expand((x+1)**2 - 1)

\text{expand}\; (x+1)^2 - 1

(12)

Factorization¶

For polynomials, factorization is the reverse of expansion. sp.factor() will convert a polynomial to a product of irreducible factors over the rational numbers. For example:

sp.factor(x**2 + 2*x + 1)

\text{factor}\; x^2 + 2x + 1

(13)

You can actually pass any expression to sp.factor() and sp.expand() and they will try to factor or expand it, even it is not a pure polynomial.

sp.factor(sp.sin(x)**2 + sp.cos(x)**2)

\text{factor}\; \sin^2 x + \cos^2 x

(14)

sp.expand(sp.sin(x)**2 + sp.cos(x)**2)

\text{expand}\; (\sin x + \cos x)^2

(15)

Collection¶

The sp.collect() function will combine terms in an expression that share a common factor, and it takes two arguments: the expression to be combined, and the factor to be combined. For example:

sp.collect(x*y + x - 3 + 2*x**2 - z*x**2 + x**3, x)

\text{collect}\; xy + x - 3 + 2x^2 - zx^2 + x^3 \; \text{for $x$}

(16)

Cancellation¶

The sp.cancel() function simplifies a rational expression into its canonical form $\frac{p}{q}$ , where $p$ and $q$ are expanded polynomials that share no common factors. In this canonical form, the leading coefficients of $p$ and $q$ are integers, meaning they have no denominators.

sp.cancel((x**2 + 2*x + 1)/(x**2 + x))

\text{cancel}\; \frac{x^2 + 2x + 1}{x^2 + x}

(17)

Partial Fraction Decomposition¶

The sp.apart() function will decompose a rational expression into a product of linear combinations of monomials; i.e., the partial fraction decomposition. For example:

sp.apart((4*x**3 + 21*x**2 + 10*x + 12)/(x**4 + 5*x**3 + 5*x**2 + 4*x))

\text{decompose}\; \frac{4x^3 + 21x^2 + 10x + 12}{x^4 + 5x^3 + 5x^2 + 4x}

(18)

Trigonometric Simplification¶

The sp.trigsimp() function can simplify trigonometric expressions using trigonometric identities.

sp.trigsimp(sp.sin(x)**2 + sp.cos(x)**2)

\text{trig simplify}\; \sin^2 x + \cos^2 x

(19)

This can also do for hyperbolic trigonometric functions:

sp.trigsimp(sp.sinh(x)**2 + sp.cosh(x)**2)

\text{trig simplify}\; \sinh^2 x + \cosh^2 x

(20)

The reverse of trigonometric simplification is sp.expand_trig(), which applies the sum or double angle identities to expand trigonometric functions.

sp.expand_trig(sp.sin(x + y))

\text{trig expand}\; \sin(x + y)

(21)

There are definitely more functions to simplify expressions. You can check the SymPy tutorial of Simplification and API of Simplification for more details.

Exercise : Simplify the following expressions.

$\displaystyle \frac{x^2 - 1}{x - 1}$
$\tan^2 x + 1$
$(x + y)^2 - (x^2 + 2xy + y^2)$
$\displaystyle \frac{x^3 - 3x^2 + 3x - 1}{(x - 1)^3}$
$\displaystyle \frac{\sinh(2x)}{2\sinh(x)\cosh(x)}$
$\displaystyle \frac{e^{2\ln(x)} - 1}{e^{\ln(x)} - 1}$

Calculus¶

Calculus is probably the first nightmare most students face when they step into college. As chemists, we might not always have the strongest grip on it—but luckily, we don’t have to! With SymPy, we can let Python handle all those tedious derivatives, integrals, and limits for us. In this section, we’ll see how symbolic calculus works in SymPy and how it can make life a lot easier when dealing with real chemical and physical problems.

Derivatives¶

Finding derivatives is a very common task in chemistry and physics. It is also easy to perform in SymPy. You can use the sp.diff() function to find the derivative of an expression. This function takes two arguments: the expression to be differentiated, and the variable to be differentiated with respect to. For example:

sp.diff(sp.sin(x), x)

\frac{d}{dx} \sin x

(22)

We can also take derivatives of higher order, by either passing the variable to be differentiated and the order of the derivative as separate arguments:

sp.diff(x**3 + 2*x + 1, x, 3)

\frac{d^3}{dx^3} (x^3 + 2x + 1)

(23)

Or by passing the variable to be differentiated for multiple times:

sp.diff(x**3 + 2*x + 1, x, x, x)

\frac{d}{dx} \frac{d}{dx} \frac{d}{dx} (x^3 + 2x + 1)

(24)

For partial derivatives, the implementation is similar—just pass each variable to be differentiated and the order of the derivative as separate arguments:

sp.diff(sp.exp(x*y*z), x, 2, y, 1, z, 3)

\frac{\partial^6}{\partial x^2 \partial y \partial z^3} e^{x y z}

(25)

For clarity, you can collect each variable and its order in a tuple:

sp.diff(e**(x*y*z), (x, 2), (y, 1), (z, 3))

diff() can also be called as a method:

expr = sp.sin(x**2) * sp.exp(x*y)
expr.diff(x, 2, y)

\frac{\partial^3}{\partial x^2 \partial y} \sin(x^2) e^{xy}

(26)

We can also make an unevaluated derivative by using the sp.Derivative class. It has the same syntax as sp.diff() (but you can’t use it as a method).

sp.Derivative(sp.sin(x**2), x, 2, y)

To evaluate an unevaluated derivative, you can use the doit() method:

deriv_expr = sp.Derivative(expr, x, 2, y)
deriv_expr.doit()

These unevaluated derivative objects are useful when you want to postpone evaluation or display the derivative symbolically. They are also used in cases where SymPy cannot directly compute the derivative, such as when dealing with differential equations that include undefined functions.

Derivatives of any order can be defined by using a tuple (x, n), where n specifies the order of differentiation with respect to x.

deriv_any_order = (a*x + b)**m
deriv_any_order.diff((x, n))

For more details, you can check the SymPy documentation of diff() and SymPy documentation of Derivatives.

Exercise : Find the following derivatives.

$\displaystyle \frac{d}{dt} \left( \frac{2t}{\sqrt{t} - e^{-t}} + \frac{3t}{\sqrt{t} + e^{-t}} \right)$
$\displaystyle \frac{\partial^2}{\partial \theta^2} \left( \left( \frac{1 + \sin\theta}{1 - \cos\phi} + \cosh r \right)^2 \right)$
$\displaystyle \frac{\partial^2}{\partial x \partial y} \left( \sqrt{x^5 + y^5} \csc(e^{x+y}) \right)$

Exercise : Given the Morse potential

V_\text{Morse}(r) = D_e \left( 1 - e^{-a \left( r - r_0 \right)} \right)^2

(27)

find the first derivative of $V_\text{Morse}$ with respect to $r$ . Then use subs() to set $r = r_0$ to verify that the derivative is zero at the equilibrium distance.

Integrals¶

There are two kinds of integrals: definite and indefinite. To compute an indefinite integral, we can use the sp.integrate() function. The first argument is the expression to be integrated, and the second argument is the variable to be integrated with respect to. For example:

sp.integrate(sp.sin(x), x)

\int \sin x \, dx

(28)

Note that the integration constant is NOT included in the result.

To compute a definite integral, just pass the variable to be integrated and the lower and upper bounds of the integration interval as the second argument, wrapped in parentheses (tuple):

sp.integrate(sp.sin(x), (x, 0, 2*sp.pi))

\int_0^{2\pi} \sin x \, dx

(29)

For integrals containing $\infty$ , you can pass sp.oo as the corresponding bound—this is $\infty$ in SymPy—a good imitation, isn’t it?

sp.integrate(sp.exp(-x), (x, 0, sp.oo))

\int_0^\infty e^{-x} \, dx

(30)

Multiple integrals can be computed by passing multiple tuples:

sp.integrate(sp.exp(-x**2 - y**2), (x, -sp.oo, sp.oo), (y, -sp.oo, sp.oo))

\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} e^{-x^2 -y^2} \, dx \, dy

(31)

If sp.integrate() is unable to compute an integral, it will return an unevaluated sp.Integral object.

sp.integrate(x**x, x)

\int x^x \, dx

(32)

To evaluate an unevaluated sp.Integral, you can use the doit() method, similar to how we did for sp.Derivative.

integral = sp.Integral(x**2 * sp.exp(-x**3), x)
integral

integral.doit()

SymPy can even output results as special functions or piecewise functions:

sp.integrate(sp.sin(x**2), x)

\int \sin x^2 \, dx

(33)

Here $S(x)$ is the Fresnel integral of sine.

Also see this:

sp.integrate(x**y*sp.exp(-x), (x, 0, sp.oo))

\int_0^\infty x^y e^{-x} \, dx

(34)

Here the output is a piecewise function, with $\mathrm{re}(y)$ the real part of $y$ .

Integrals can also be evaluated numerically. This can be done by using the method evalf() on the sp.Integral object.

integral2 = sp.Integral(x * sp.exp(x**2), x)
integral2.evalf()

\int_0^1 x e^{x^2} \, dx

(35)

One thing to be aware of is that numerical integration in SymPy can be quite slow. The evalf() and subs() methods can provide accurate results with any desired level of precision, but if you value speed over accuracy, especially when working with large datasets, it’s better to use another library, such as NumPy or SciPy, to handle these calculations. One way to do this is by using sp.lambdify() to convert a SymPy expression into a function that can be evaluated efficiently by numerical libraries, as introduced earlier. For more methods, check the SymPy documentation on Numerical Computation.

Exercise : Evaluate the following integrals.

$\displaystyle \int_0^1 \frac{1 - \cos 2x}{2x^3} \, dx$
$\displaystyle \int \frac{1}{x^5 + 1} \, dx$
The so-called elliptic integral: $\displaystyle \int \frac{d\theta}{ \sqrt{1 - k^2 \sin^2 \theta} }$

Exercise : Integral with variable limits can also be computed in SymPy. Verify the following expression, which is called the general Leibniz integral rule:

\frac{d}{dx} \left( \int_{a(x)}^{b(x)} f(x, t) \, dt \right) = f(x, b(x)) \frac{d}{dx} b(x) - f(x, a(x)) \frac{d}{dx} a(x) + \int_{a(x)}^{b(x)} \frac{\partial}{\partial x} f(x, t) \, dt.

(36)

You should define your function $f(x, t)$ , $a(x)$ , and $b(x)$ as follows:

f = sp.Function('f')(x, t)
a = sp.Function('a')(x)
b = sp.Function('b')(x)

Limits¶

Limits can be computed using the sp.limit() function. The first argument is the expression to be limited, and the second argument is the variable to be limited with respect to, and the third argument is the limit point.

sp.limit(sp.sin(x)/x, x, 0)

\lim_{x \to 0} \frac{\sin x}{x}

(37)

Use sp.limit() instead of subs() when evaluating expressions at singular points. Although SymPy provides objects sp.oo to represent infinity ( $\infty$ ), direct evaluation with them isn’t always reliable because it doesn’t account for factors like rates of growth. For example, operations such as $\infty - \infty$ or $\frac{\infty}{\infty}$ return $\mathrm{NaN}$ (not-a-number, nan). You can compare these two results below:

expr_lim = x**2/sp.exp(x)
expr_lim.subs(x, sp.oo)

\lim_{x \to \infty} \frac{x^2}{e^x}

(38)

sp.limit(expr_lim, x, sp.oo)

Similar to sp.Derivative and sp.Integral, the function sp.limit() also has an unevaluated form called sp.Limit. You can evaluate it explicitly by calling doit().

lim = sp.Limit((sp.cos(x) - 1)/x, x, 0)

lim.doit()

To evaluate a one-sided limit, provide '+' or '-' as the fourth argument in sp.limit().

sp.limit(1/x, x, 0, '+')

\lim_{x \to 0^+} \frac{1}{x}

(39)

From the other side:

sp.limit(1/x, x, 0, '-')

\lim_{x \to 0^-} \frac{1}{x}

(40)

Exercise : In relativity, the observed length of a moving object, such as a rocket, depends on its velocity relative to the observer—this is called the length contraction or Lorentz contraction.

If the rocket’s length at rest is $L_0$ , then at speed $v$ its length appears as:

L = L_0 \sqrt{1 - \frac{v^2}{c^2}}.

(41)

Here, $c$ represents the speed of light in a vacuum.

How does $L$ change as $v$ increases? Find the limit $\displaystyle \lim_{v \to c^-} L$ . Why is the left-hand limit required in this case?

Question is adapted from Hass, J. R.; Heil, C. E.; Weir, M. D. Thomas’ Calculus, 14th ed.; Pearson: New York, 2017. ISBN 9780134438986.

Series Expansion¶

Power Series¶

Power series are really useful for approximating functions. They are defined as:

\sum_{n=0}^{\infty} a_n x^n = a_0 + a_1 x + a_2 x^2 + \cdots

(42)

For well-behaved (analytic) functions, we can use power series to precisely represent them. You can also perform a power series expansion to get arbitrary-precision approximations.

In SymPy, power series expansion can be computed using the sp.series(expr, x=None, x0=0, n=6, dir='+') function. The first argument is the expression to be expanded, and the second argument is the variable to be expanded with respect to. The third argument is the expansion point, and the fourth argument is the maximum degree of the expansion. The fifth argument is the direction of the expansion, which can be '+' or '-'. If an infinite expansion is desired, set the third argument to sp.oo. For negative infinity, set dir=-`.

sp.series(sp.exp(x), x)

Here, $O(x^n)$ denotes the asymptotic order of the remainder, meaning that as $x \to x_0$ , the error term grows no faster than a constant multiple of $(x - x_0)^n$ . If you don’t want to see the remainder, you can use the method removeO() to remove it (O is letter O but not number 0!).

sp.series(sp.exp(x), x).removeO()

For more details, check the SymPy documentation of series() (it’s not well-written, though...).

Fourier Series¶

Fourier series provide a way to represent a periodic function as a sum of trigonometric functions. There are several conventions for normalizing Fourier series; SymPy adopts the following form:

\frac{a_0}{2} + \sum_{n=1}^{\infty} \left[ a_n \cos \left( \frac{2\pi n x}{L} \right) + b_n \sin \left( \frac{2\pi n x}{L} \right) \right]

(43)

This represents the Fourier series of a function $f(x)$ over the interval $(a, b)$ , where $L = b - a$ . The coefficients $a_n$ and $b_n$ are defined as:

a_0 = \frac{2}{L} \int_a^b f(x)\,dx

(44)

a_n = \frac{2}{L} \int_a^b f(x) \cos \left( \frac{2\pi n x}{L} \right)\,dx

(45)

b_n = \frac{2}{L} \int_a^b f(x) \sin \left( \frac{2\pi n x}{L} \right)\,dx

(46)

It turns out that it is not strictly necessary for $f(x)$ to be periodic, as the Fourier series can still be defined and analyzed over a finite interval.

To compute a Fourier series expansion in SymPy, use the sp.fourier_series(f, limits=None) function. The first argument is the expression to be expanded, and the second argument is a tuple specifying the variable along with its lower and upper limits for the expansion interval. If no limits are provided, the default interval is set to $(-\pi, \pi)$ .

For more details, check the SymPy documentation of fourier_series().

Exercise : Find the 10th-degree power series expansion of the following functions using SymPy:

$\displaystyle \frac{1}{1 + e^{-x}}$
$\displaystyle \sin \left(\frac{1}{x^2+1} \right) \ \text{at}\ x_0 = -1$
$\displaystyle \int \cos x^2 \, dx$

Exercise : Find the Fourier series expansion of the following square wave defined over one period $(0, 2)$ :

p = sp.Piecewise(
    (1, (x >= 0) & (x < 1)),
    (0, (x >= 1) & (x < 2))
)

Solving Equations¶

Probably the most important thing we want to do with a computer is to solve equations. SymPy provides a variety of functions for solving equations. We will briefly cover the process of solving algebraic equations and leave differential equations and other types of equations in Lesson A2.

A thing to be noticed is that the equation in SymPy is not represented by the assignment operator (=), nor the equality operator (==). Instead, the equation needs to be created as an object of the sp.Eq class. For example, you can create the equation

x^2 + 2x + 1 = 0

(47)

by using the following code:

sp.Eq(x**2 + 2*x + 1, 0)

There are two kinds of equation solvers in SymPy: sp.solve() and sp.solveset(). The main difference is that sp.solve() will provide exact forms of symbolic representations of the solutions, while sp.solveset() will provide a set of solutions. They sound similar, but in practice, the result of sp.solve() can use the method subs() but sp.solveset() cannot; also, sp.solveset() can provide more than one solution that sp.solve() sometimes cannot, especially for periodic solutions. Besides these, their syntax is basically the same.

Let’s see some examples:

sp.solve(eq_1, x)

sp.solveset(eq_1, x)

sp.solve(sp.Eq(sp.sin(x), 0), x)

sp.solveset(sp.Eq(sp.sin(x), 0), x)

See? For simple equations, sp.solve() and sp.solveset() are the same, but for periodic solutions, sp.solveset() can catch all solutions, while sp.solve() will only give solutions in one period (but there are some ways that you can generate a full set of solutions from this one period, which we will not cover in this lesson).

For obtaining solutions in a specific interval, you can use the sp.solveset() function with the sp.Interval(start,end) class. For example,

sp.solveset(sp.Eq(sp.sin(x), 0), x, sp.Interval(-1, 1))

Note that the interval is closed, so the result will include the endpoints. You can use left_open and right_open arguments to specify whether the interval is open on the left or right side.

sp.solveset(sp.Eq(sp.sin(x), 0), x, sp.Interval(0, 3*sp.pi, left_open=True, right_open=True))

You can also designate a specific set for the solution by using the sp.S class. Some useful sets are shown in the table below.

Set	Description
`sp.S.Rationals`	Rational numbers
`sp.S.Reals`	Real numbers
`sp.S.Complexes`	Complex numbers
`sp.S.Integers`	Integers
`sp.S.Naturals`	Natural numbers
`sp.S.Naturals0`	Non-negative integers
`sp.S.Positive`	Positive numbers
`sp.S.Negative`	Negative numbers

For example, you can solve the equation $x^3 = 1$ for $x$ in the set of real and complex numbers by using the following code:

sp.solveset(sp.Eq(x**3, 1), x, sp.S.Reals)

sp.solveset(x**3, x, sp.S.Complexes)

A shortcut here is, you can pass a simple expression to sp.solve() and sp.solveset(), without creating an sp.Eq object. In this case, solvers will automatically interpret the expression as an equation, with the right-hand side equal to 0.

sp.solve(x**2, x)

sp.solveset(x**2, x)

End-of-Lesson Problems¶

Problem 1: The Most Probable Speed of Gaseous Molecules¶

Remember that in the introduction section, we mentioned that the most probable speed $v_\text{mp}$ of gaseous molecules can be interpreted from the Maxwell-Boltzmann distribution

f(v) = 4\pi \left( \frac{m}{2\pi k_B T} \right)^\frac{3}{2} v^2 e^{-\frac{mv^2}{2k_BT}}

(48)

by finding the maximum of $f(v)$ . Use methods learned in this lesson to find the expression of $v_\text{mp}$ .

If you find any weird behavior of your result (e.g., output is not exact; multiple solutions; etc.), please redefine the variables using sp.symbols() with appropriate assumptions and try again.

Also, we can find the mean speed $\bar{v}$ of gaseous molecules. The mean speed is simply the expected value of the speed distribution, which can be computed by integrating the distribution multiplying $v$ over the range of speeds:

\bar{v} = \int_0^\infty v f(v)\, dv

(49)

Find the expression of $\bar{v}$ .

Problem 2: Archimedes’ Principle¶

In a static, homogeneous fluid of density $\rho$ , with gravity $g$ acting downward, the hydrostatic pressure at point $\mathbf{r} = (x,y,z)$ is

p = p_0 - \rho g z

(50)

where $p_0$ is the pressure at $z=0$ and $z$ is the vertical coordinate (positive upward). A solid sphere of radius $R$ is centered at the origin and is completely submerged in the fluid. Since our system is spherically symmetric, using a spherical coordinate system can significantly simplify our work. Find the expression of the pressure at the surface of the sphere using the coordinate transformation formula

(x,y,z) = \left( R\sin{\theta}\cos{\phi},\; R\sin{\theta}\sin{\phi},\; R\cos{\theta} \right)

(51)

where $\theta \in [0, \pi]$ and $\phi \in [0, 2\pi)$ .

Using the outward unit normal $\mathbf{n}=\dfrac{\mathbf{r}}{R}$ and the surface element $dS = R^2\sin{\theta}\,d\theta\,d\phi$ , we can find the force acting on the surface element by the definition $d\mathbf{F} = -p\,\mathbf{n}\,dS$ , and then integrate it on the whole surface to get the total force:

\mathbf{F} = -\iint_{S} p\,\mathbf{n}\,dS = -R^2 \int_0^{2\pi} \int_0^\pi p\,\mathbf{n}\,\sin{\theta}\,d\theta\,d\phi

(52)

We haven’t covered how to define and use vectors in SymPy yet. As an alternative, we can integrate each of its components separately.

Compare the result with $\mathrm{F} = \rho g V$ where $V$ is the volume of the sphere, which is the expression derived from Archimedes’ principle.

SymPy Documentation
Charles J. Weiss’s Scientific Computing for Chemists with Python
GenAI for making paragraphs and codes(・ω< )★
And so many resources on Reddit, StackExchange, etc.!