Object-oriented programming, symbolic computation

Object-oriented programming

In programming there are various approaches; the Python documentation (docs.python.org) mentions, for example:

1. procedural: a list of instructions on what to do (e.g.: C, Pascal)

2. declarative: we describe what we want, and the programming language carries it out (e.g.: SQL)

3. functional: programming based on functions (e.g.: Haskell )

4. object-oriented: the program is based on objects that have properties, functions, … (e.g.: Java, Smalltalk)

Python is an object-oriented programming language, but it does not force us to use objects in every case. As we will see later, object-oriented programming has many advantages, but it can often be cumbersome and would make a program unnecessarily complex. For this reason, we avoid explicit object-oriented programming whenever possible.

Object-oriented programming in Python is based on classes, and objects are instances of a class. We will look at only some of the basics of object-oriented programming (so that you can more easily understand the code of other authors and adapt it to your own needs).

A class is defined with the class statement (documentation):

class ClassName:
    '''docstring'''
    [expression 1]
    [expression 2]
    .
    .

where the class name (i.e., ClassName) is, per PEP8, written with a capital first letter. If the name consists of several words, each one is capitalized (the so-called CamelCase convention).

Let’s look at an example:

class Rectangle:
    """Class for a rectangle object"""

    def __init__(self, width=1, height=1): # this is the object constructor. It runs when we call Rectangle(width=1, height=4)
        self.width = width
        self.height = height # height is an attribute of the object

    def area(self):
        return self.width * self.height

    def set_width(self, width=1):
        self.width = width

Before we go into the details of understanding the code, let’s create an instance of the class (i.e., an object):

my_rectangle = Rectangle()

Functions defined inside a class are called methods when we call them on objects.

In the example above, we use the area method like this:

my_rectangle.area()

1

The name self is a reference to the instance of the class (the object that will be created). Names defined inside the class become attributes of the object.

An example of the height attribute:

my_rectangle.height

1

Where does the result 1 come from? When we create an object, the initialization function __init__() runs first; its arguments are the arguments we pass to the class.

Example:

your_rectangle = Rectangle(height=5, width=3)
your_rectangle.area()

15

We have also prepared a method that changes the width attribute:

my_rectangle.set_width(width=100)
my_rectangle.area()

100

We can also change attributes directly, but we usually avoid this (because of the possibility of errors and misuse).

Example:

my_rectangle.height = 3
my_rectangle.area()

300

Inheritance

An important property of classes is inheritance; the name itself says it all: just as we humans inherit from our parents, the same holds for classes. Each class (class) can thus inherit the properties of some other class (documentation). It can even have several parents (we will not go into these details here).

The syntax of a class that inherits is:

class Child(Parent):
    [expression]
    .
    .

Note: even if we do not define a parent for a class, the class inherits from object.

An example where the new class Square inherits from an existing one (Rectangle):

class Square(Rectangle):
    "Square class"

    def __init__(self, width=1):
        # call the initialization of the Rectangle class
        super().__init__(width=width, height=width)

    def set_width(self, width):
        self.width = width
        self.height = width

Now let’s look at how it is used:

my_square = Square(width=4)

The Square class does not define a method for computing the area, but it inherited it from the Rectangle class, and therefore has a method for computing the area:

my_square.area()

16

If we change the width, the area changes accordingly:

my_square.set_width(5)
my_square.area()

25

Example of inheriting from the list class

First we prepare a list:

my_list = list([1,2,3])
my_list

[1, 2, 3]

If we want to add a value to the list, we use the append method (this is a method that objects of type list have):

my_list.append(1)

Then we plot the list (first we import matplotlib):

import matplotlib.pyplot as plt
plt.plot(my_list)
plt.show()

/tmp/ipykernel_3440/2547190459.py:3: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  plt.show()

If we do something like this often, then it is better to prepare our own List class that inherits from list and add a method draw for plotting:

class List(list):

    def draw(self):
        plt.plot(self, 'r.', label='Long text')
        plt.legend()
        plt.ylim(-5, 5)
        plt.show()

An instance of the object is:

my_list_obj = List([1, 2, 3])
type(my_list_obj)

__main__.List

my_list_obj

[1, 2, 3]

Although we did not define the append method, the new class inherited it from the list class:

my_list_obj.append(1)
my_list_obj

[1, 2, 3, 1]

It also has a method for plotting:

my_list_obj.draw()

/tmp/ipykernel_3440/705516469.py:7: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  plt.show()

Symbolic computation with SymPy

The term symbolic computation means that we solve mathematical expressions by machine in the form of abstract symbols (and not numerically). Machine symbolic computation helps us when we are interested in the result in symbolic form and the expressions are too large to solve by hand on paper. We also resort to machine solving to reduce the chance of error (in extensive computations people can make mistakes).

Symbolic computation is by no means a substitute for numerical methods!

We will look at some of the basics; some recommended additional resources are:

• J.R. Johansson Scientific python lectures,

• SymPy - the official documentation of the module,

• An excellent article by some of the SymPy authors: Meurer et. al, 2017.

SymPy is one of the computer algebra systems (CAS), which, in addition to its capabilities, has the advantage of being written entirely in Python (the alternative Python package Sage, for example, is not written entirely in Python).

Some dedicated commercial tools:

• Mathematica

• Maple

First we import the SymPy module; typically the package is imported as sym:

import sympy as sym

You may notice that SymPy is also imported with from sympy import *. As a rule, we avoid this, since it unnecessarily fills the program’s base namespace with SymPy names. In the latter case, the package’s functions (e.g., sympy.Sum) are accessed directly (e.g., Sum), which can be appealing but starts to bother us when we add other packages (e.g., numpy), which, in addition to confusion, can lead to functions with the same name being “overwritten”.

In order to get nicely formatted LaTeX output, we use:

sym.init_printing()

Defining variables and numerical computation

We define variables like this:

x, y, k = sym.symbols('x, y, k')

We can check the type:

type(x)

sympy.core.symbol.Symbol

We notice that the variable x is now a Symbol from the sympy package.

Now we can make a simple computation:

x**y + k

We can also name a function; here is an example where we use the sine function and the constant \(\pi\):

f = sym.sin(1.2*sym.pi + x)**2
f

The reader may be wondering why we need a new function sympy.sin(), since we already have the one from the numpy package! The reason is that symbolic computation requires a completely different treatment than numerical computation, and therefore the code behind it is entirely different.

If we want to write an equation, i.e., to equate one expression to another, we do it like this:

equation = sym.Eq(sym.sin(k*x),0.5)
equation

Undefined mathematical functions are written as (documentation):

g = sym.Function('g')

Now we can, for example, define a differential equation:

sym.Eq(g(x).diff(x), x)

When defining variables, we can add assumptions:

x = sym.Symbol('x', positive=True)

x.is_positive

True

The assumptions are then taken into account in the computation. In general, we know that \(\sqrt{(x^2)}\ne x\), but if \(x\) is positive, then \(\sqrt{(x^2)}= x\), and sympy, given the assumptions, computes the correct result:

sym.sqrt(x**2)

With the assumptions0 method we look at an object’s assumptions:

x.assumptions0

{'commutative': True,
 'complex': True,
 'extended_negative': False,
 'extended_nonnegative': True,
 'extended_nonpositive': False,
 'extended_nonzero': True,
 'extended_positive': True,
 'extended_real': True,
 'finite': True,
 'hermitian': True,
 'imaginary': False,
 'infinite': False,
 'negative': False,
 'nonnegative': True,
 'nonpositive': False,
 'nonzero': True,
 'positive': True,
 'real': True,
 'zero': False}

SymPy knows the types of numbers (documentation):

• Real real numbers,

• Rational rational numbers,

• Complex complex numbers,

• Integer integers.

These types are important, as they can affect the way of solving and the solution.

Rational numbers

Above we have already defined real numbers. Let’s now look at rational numbers with an example:

r1 = sym.Rational(4, 5)
r2 = sym.Rational(5, 4)

r1

Some computations:

r1+r2

r1/r2

Complex numbers

The imaginary number is written with I (this is different from numpy, where the imaginary part is defined with j):

1+1*sym.I

sym.I**2

Numerical computation

In symbolic computation we first solve analytical expressions, simplify them, etc., and then we often also want to compute a concrete result.

Let’s look at an example:

x = sym.symbols('x')
f = sym.exp(x**x**x + sym.pi)
f

If we now want to use a value instead of \(x\), e.g., 0.5, we do it with the subs() method (documentation):

f.subs(x, 0.5)

Above we used the constant \(\pi\) (documentation); some commonly used constants are:

• sympy.pi for the number \(\pi\),

• sympy.E for the natural number \(e\),

• sympy.oo for infinity.

As we saw above, subs() performs the substitution and then the simplifications that are obvious; it did not compute the number \(\pi\) in rational form. We must explicitly request this with the method:

• evalf (evaluate function, documentation) or

• N,

both of which have the argument n (number of digits).

Let’s look at an example:

f.subs(x, 2).evalf(n=80)

It is similar with N:

sym.N(f.subs(x, 2), n=80)

In passing, we have shown that, provided the expression contains no floating-point numbers, we can display the result to arbitrary precision (documentation).

In the subs function we can also use a dictionary. Example:

x, y = sym.symbols('x, y')
parameters = {x: 4, y: 10}

(x + y).subs(parameters)

or a list of tuples:

(x + y).subs([(x, 4), (y, 10)])

Similarly, the sympy.evalf() method has a subs argument that accepts a dictionary of substitutions, for example:

(y**x).evalf(subs=parameters)

The sympy.subs method can also replace a symbol (or expression) with another expression:

(x + y).subs(x, y + sym.oo)

SymPy and NumPy

We often connect sympy with numpy. As an example, let’s look at how we would compute the expression:

x = sym.symbols('x')
f = sym.sin(x) + sym.exp(x)
f

efficiently and numerically at a thousand values of \(x\).

First we import the numpy package:

import numpy as np

We prepare a numerical array of values:

x_vec = np.linspace(0, 10, 1000, endpoint=False)
x_vec[:10]

array([0.  , 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09])

Based on what is written above and our prior knowledge, we use a list comprehension:

y_vec = np.array([f.evalf(subs={x: value}) for value in x_vec])

We notice that this is slow, so let’s measure the time required:

%%timeit -n1
y_vec = np.array([f.evalf(subs={x: value}) for value in x_vec])

78.1 ms ± 593 μs per loop (mean ± std. dev. of 7 runs, 1 loop each)

Using the lambdify function

A significantly faster approach is to use the lambdify approach, where a compiled function is prepared, optimized for numerical execution. The syntax of the sympy.lambdify() function is (documentation):

sympy.lambdify(symbols, function, modules=None)

where the arguments are:

• symbols the symbols used in function that are replaced with numerical values,

• function represents a sympy function,

• modules represents which package the compiled form is prepared for. If numpy is installed, it is the default for that module.

Example of use:

f_fast = sym.lambdify(x, f, modules='numpy')
y_vec_fast = f_fast(x_vec)

Let’s check the speed:

%%timeit -n100
f_fast(x_vec)

19.5 μs ± 1.56 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)

We observe roughly a 10,000-fold speedup!

Let’s also look at an example of using a function of several variables:

f_fast2 = sym.lambdify((x, y), (x + y + sym.pi)**2, 'numpy')
x = np.linspace(0, 10, 10)
y = x
f_fast2(x, y)

array([  9.8696044 ,  28.77051002,  57.54795885,  96.20195089,
       144.73248614, 203.1395646 , 271.42318627, 349.58335115,
       437.62005924, 535.53331054])

Graphical display

SymPy has data plotting based on matplotlib. The plotting is more limited compared to matplotlib and we use it for simple plots. (documentation).

We will look at simple examples related to the sympy.plotting.plot function; first we import the function:

from sympy.plotting import plot

The syntax for using the sympy.plotting.plot() function (documentation) is:

plot(expression, range, **kwargs)

where the arguments are:

• expression is a mathematical expression or several expressions,

• range is the range of the plot (the default range is (-10, 10)),

• **kwargs are keyword arguments, i.e., a dictionary of various options.

The function returns an instance of the sympy.Plot() object.

A minimal example of a single function:

x = sym.symbols('x')
plot(x**2)

/opt/hostedtoolcache/Python/3.11.15/x64/lib/python3.11/site-packages/sympy/plotting/backends/matplotlibbackend/matplotlib.py:309: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
  self.plt.show()

<sympy.plotting.backends.matplotlibbackend.matplotlib.MatplotlibBackend at 0x7fdc1ca0b110>

Note: the last line points to the result in the form of a Plot instance; we hide the output of the object <sympy.plotting.plot.Plot at 0x...> by using a semicolon:

plot(x**2);

but the figure is still displayed.

Some common arguments are:

• show displays the figure (default True),

• line_color the color of the plot,

• xscale and yscale the display mode (options: linear or log),

• xlim and ylim the display limits for the axes (a tuple of two (min, max) values).

Let’s prepare two figures where the \(y\) axis is logarithmic and the plot range is from 1 to 5:

plot1 = plot(x**2, (x, 1, 5), show=False, legend=True, yscale='log')
plot2 = plot(30*sym.log(x), (x, 1, 5), show=False, line_color='C2', legend=True, yscale='log')

Now we extend the first figure with the second and display the result:

plot1.extend(plot2)
plot1.show()

Parametric plot

Similarly, we use the function sympy.plotting.plot_parametric for a parametric plot (documentation):

plot_parametric(expr_x, expr_y, range, **kwargs)

where the new arguments are:

• expr_x and expr_y the definitions of the position of coordinates \(x\) and \(y\),

• **kwargs a dictionary of options.

We import the function:

from sympy.plotting import plot_parametric

Let’s show its use with an example:

plot_parametric(sym.sin(x), sym.sin(2*x), (x, 0, 2*sym.pi));

Plotting in space

The function sympy.plotting.plot3D (documentation) for plotting in space has the syntax:

plot3d(expression, range_x, range_y, **kwargs)

where the arguments are:

• expression the definition of the surface,

• range_x and range_y the range of coordinate x and y,

• **kwargs a dictionary of options.

We import the function:

from sympy.plotting import plot3d

Let’s show its use with an example:

x, y = sym.symbols('x y')
plot3d(x**2 + y**2, (x, -5, 5), (y, -5, 5));

For other plots, see the documentation.

Algebra

In this section we will look at some of the basics of using SymPy for algebraic operations.

Using expand and factor

Let’s define a mathematical expression:

x = sym.symbols('x')
f = (x+1)*(x+2)*(x+3)
f

and now expand it (see the documentation):

aa = sym.expand(f)
aa

If we want to look at the coefficients in front of x, we do it with the coeff() method:

aa.coeff(x)

The arguments of the function define what kind of expansion we want (documentation). If we want, e.g., a trigonometric expansion, then we use trig=True:

a, b = sym.symbols('a, b')
sym.expand(sym.sin(a+b))

sym.expand(sym.sin(a+b), trig=True)

The inverse operation of expansion is factorization (documentation):

sym.factor(x**3 + 6 * x**2 + 11*x + 6)

If we are interested in the individual factors, then we do this with the sympy.factor_list function:

sym.factor_list(x**3 + 6 * x**2 + 11*x + 6)

Simplifying expressions with simplify

The sympy.simplify() function (documentation) tries to simplify expressions into simpler ones (e.g., by canceling variables).

For special purposes we can also simplify with:

• sympy.trigsimp,

• sympy.powsimp,

• sympy.logcombine.

For more, see the documentation.

Examples of simplification:

sym.simplify((x+1)*(x+1)*(x+3))

sym.simplify(sym.sin(a)**2 + sym.cos(a)**2)

sym.simplify(sym.cos(x)/sym.sin(x))

Using apart and together

We use these functions to work with fractions:

f1 = 1/((1 + x) * (5 + x))
f1

We perform partial fraction decomposition with the sympy.apart() function (documentation):

f2 = sym.apart(f1, x)
f2

and then back again in the opposite direction with the sympy.together() function:

sym.together(f2)

In the latter case we arrive at a similar result with sympy.simplify():

sym.simplify(f2)

Differentiation

Differentiation is in principle a relatively simple mathematical operation, which we perform with the sympy.diff() function (documentation):

Let’s prepare an example:

x, y, z = sym.symbols('x, y, z')
f = sym.sin(x*y) + sym.cos(y*z)
f

Let’s differentiate it with respect to \(x\):

sym.diff(f, x)

or also

f.diff(x)

Higher-order derivatives are defined like this:

sym.diff(f, x, x, x)

or (the same result a bit differently):

sym.diff(f, x, 3)

The derivative with respect to several variables \(\frac{d^3f}{dx\,dy^2}\) is performed like this:

sym.diff(f, x, 1, y, 2)

Integration

The integrate function can be used for indefinite integration (documentation):

integrate(f, x)

or for definite integration:

integrate(f, (x, a, b))

where the arguments are:

• f the function we integrate,

• x the variable with respect to which we integrate,

• a and b the limits of integration.

An example of indefinite integration:

x = sym.symbols('x')
f = sym.sin(x*y) + sym.cos(y*z)
sym.integrate(f, x)

\[\begin{split}\displaystyle x \cos{\left(y z \right)} + \begin{cases} - \frac{\cos{\left(x y \right)}}{y} & \text{for}\: y \neq 0 \\0 & \text{otherwise} \end{cases}\end{split}\]

We notice that sympy correctly takes into account the possibility that \(y=0\).

Another example of definite integration:

sym.integrate(f, (x, -1, 1))

An example where the limits are at infinity (we use the constant for infinity sympy.oo):

sym.integrate(sym.exp(-x**2), (x, -sym.oo, sym.oo))

Sum and product of a series

We define the sum of a series with the sympy.Sum() function (documentation):

sympy.Sum(expression, (var, start, end))

where the arguments are:

• expression the expression we are summing,

• var, start and end the variable that increases from start to end (end is included).

An example of the sum of a series:

n = sym.Symbol('n')
f = sym.Sum(1/x**n, (n, 1, sym.oo))
f

Only when we use the doit() method is the computation carried out:

f.doit()

\[\begin{split}\displaystyle \begin{cases} \frac{1}{x \left(1 - \frac{1}{x}\right)} & \text{for}\: \frac{1}{\left|{x}\right|} < 1 \\\sum_{n=1}^{\infty} x^{- n} & \text{otherwise} \end{cases}\end{split}\]

Let’s also look at the numerical result:

f.subs({x: 5}).evalf()

We define the product of a series similarly with the sympy.Product function (documentation):

sympy.Product(expression, (var, start, end))

where the arguments are:

• expression the expression we are multiplying,

• var, start and end the variable that increases from start to end (end is included).

Example:

f = sym.Product(1/n, (n, 1, 5))
f

f.doit()

Limits

We compute limits with the sympy.limit() function (documentation):

sympy.limit(f, x, x0)

where the arguments are:

• f the expression whose limit we are seeking,

• x the variable that approaches x0,

• x0 the limit point.

Example:

x = sym.symbols('x')
f = sym.sin(x)/x
f

sym.limit(f, x, 0)

As an example, let’s look at the use of the limit in the definition of the derivative:

\[ \frac{\mathrm{d}f}{\mathrm{d}x} = \lim_{h\rightarrow 0}\frac{f(x+h,y)-f(x,y)}{h}. \]

Let’s prepare the function f and its derivative:

x, y, z, h = sym.symbols('x, y, z, h')
f = sym.sin(x*y) + sym.cos(y*z)

The derivative of the function is:

sym.diff(f, x)

We compute the same result also using the limit:

sym.limit((f.subs(x, x+h) - f)/h, h, 0)

Taylor series

We compute Taylor series with the sympy.series() function (documentation):

sympy.series(expression, x=None, x0=0, n=6, dir='+')

where the arguments are:

• expression the expression whose series we are determining,

• x the independent variable,

• x0 the value around which we determine the series (default 0),

• n the order of the series (default 6),

• dir the direction of the series expansion (+ or -).

Example:

x = sym.symbols('x')
sym.series(sym.exp(x), x) # default values x0=0, and n=6

If we want to define a different expansion point (x0=2) and use more terms (n=8), we do it like this:

s1 = sym.series(sym.exp(x), x, x0=2, n=8)
s1

The result also includes the order of validity; in this way we can control the validity of the computation (\(\mathcal{O}\)).

Example:

s1 = sym.cos(x).series(x, 0, 5)
s1

s2 = sym.sin(x).series(x, 0, 2)
s2

The computations s1 and s2 have different orders of validity, so the product:

s1 * s2

is accurate only up to order \(\mathcal{O}(x^2)\), which sympy handles appropriately:

s3 = sym.simplify(s1 * s2)
s3

The information about the order of validity can be removed:

s3.removeO()

Linear algebra

Matrices and vectors

We define matrices and vectors with the Matrix function. While with numpy.array it is not necessary to strictly adhere to the mathematical notation of vectors and matrices, with sympy it is required.

Let’s look at an example; first we prepare the variables:

m11, m12, m21, m22 = sym.symbols('m11, m12, m21, m22')
b1, b2 = sym.symbols('b1, b2')

Then the matrix and the column vector:

A = sym.Matrix([[m11, m12],[m21, m22]])
A

\[\begin{split}\displaystyle \left[\begin{matrix}m_{11} & m_{12}\\m_{21} & m_{22}\end{matrix}\right]\end{split}\]

b = sym.Matrix([[b1], [b2]])
b

\[\begin{split}\displaystyle \left[\begin{matrix}b_{1}\\b_{2}\end{matrix}\right]\end{split}\]

Now let’s look at some typical operations; first the multiplication of a matrix and a vector:

A * b

\[\begin{split}\displaystyle \left[\begin{matrix}b_{1} m_{11} + b_{2} m_{12}\\b_{1} m_{21} + b_{2} m_{22}\end{matrix}\right]\end{split}\]

Then the dot product of two vectors (we must be careful to transpose one of the vectors):

b.T*b

\[\displaystyle \left[\begin{matrix}b_{1}^{2} + b_{2}^{2}\end{matrix}\right]\]

The determinant and the inverse matrix:

A.det()

A.inv()

\[\begin{split}\displaystyle \left[\begin{matrix}\frac{m_{22}}{m_{11} m_{22} - m_{12} m_{21}} & - \frac{m_{12}}{m_{11} m_{22} - m_{12} m_{21}}\\- \frac{m_{21}}{m_{11} m_{22} - m_{12} m_{21}} & \frac{m_{11}}{m_{11} m_{22} - m_{12} m_{21}}\end{matrix}\right]\end{split}\]

Multiplication and power of a matrix:

A*A

\[\begin{split}\displaystyle \left[\begin{matrix}m_{11}^{2} + m_{12} m_{21} & m_{11} m_{12} + m_{12} m_{22}\\m_{11} m_{21} + m_{21} m_{22} & m_{12} m_{21} + m_{22}^{2}\end{matrix}\right]\end{split}\]

A**2

\[\begin{split}\displaystyle \left[\begin{matrix}m_{11}^{2} + m_{12} m_{21} & m_{11} m_{12} + m_{12} m_{22}\\m_{11} m_{21} + m_{21} m_{22} & m_{12} m_{21} + m_{22}^{2}\end{matrix}\right]\end{split}\]

Root finding (solving equations)

We solve equations and systems of equations with the sympy.solve() function (documentation). Solving the following equations is supported:

• polynomial equations,

• transcendental equations,

• piecewise-defined equations as a combination of the above two types,

• systems of linear and polynomial equations,

• systems of equations with inequalities.

The syntax is:

sympy.solve(f, *symbols, **flags)

where the arguments are:

• f an expression or a list of expressions,

• *symbols the symbol or list of symbols we want to determine,

• **flags a dictionary of options.

Let’s look at an example:

x = sym.symbols('x')
f = sym.sin(x)
eq = sym.Eq(f, 1/2)
eq

sym.solve(eq, x)

Let’s plot the solution (notice that we found only two of the infinitely many solutions):

p1 = sym.plotting.plot(sym.sin(x), (x, -2*sym.pi, 2*sym.pi), line_color='C0', show=False, legend=True)
p2 = sym.plotting.plot(0.5, (x, -2*sym.pi, 2*sym.pi), line_color='C1', show=False, legend=True)
p1.extend(p2)
p1.show()

A quadratic equation:

a, b, c, x = sym.symbols('a, b, c, x')
sym.solve(a*x**2 + b*x + c, x)

A system of equations:

x, y = sym.symbols('x y')
sym.solve([x + y - 1, x - y - 1], [x, y])

For nonlinear systems we can also use numerical solving with the sympy.nsolve() function (documentation):

sympy.nsolve(f, [args,] x0, modules=['mpmath'], **kwargs)

where the arguments are:

• f the equation or system of equations we are solving,

• args the variables (optional),

• x0 the initial estimate (a scalar or a vector),

• modules the package used to compute the numerical value (the same logic as with the lambdify function; the default is the mpmath package),

• **kwargs a dictionary of options.

Let’s look at the example from above:

x = sym.symbols('x')
eq = sym.Eq(sym.sin(x), 0.5)
sol = sym.nsolve(eq, x, 3)
sol

Solving differential equations

Undefined functions

Before we look at differential equations, we must get to know undefined functions. We define such functions with sympy.Function documentation.

As an example, let’s look at how we define the equation \(\ddot x(t)=g\). First let’s define the undefined function:

x = sym.Function('x')

Then the symbols:

t, g = sym.symbols('t, g')

And the equation:

sym.Eq(x(t).diff(t,2), g)

Differential equations

We solve differential equations and systems of differential equations with the sympy.dsolve() function (documentation):

sympy.dsolve(eq, func=None, hint='default', simplify=True,
             ics=None, xi=None, eta=None, x0=0, n=6, **kwargs)

where the selected arguments are:

• eq the differential equation or system of differential equations,

• func the solution we are seeking,

• ics the initial and boundary conditions of the differential equation.

Let’s look at the example of a mass \(m\) sliding on a surface with friction coefficient \(\mu\); the initial velocity is \(v_0\), and the displacement \(x_0=0\).

Let’s define the symbols:

x = sym.Function('x')
t, m, mu, g, v0, x0 = sym.symbols('t m mu g v0 x0', real=True, positive=True)

Let’s define the differential equation:

eq = sym.Eq(m*x(t).diff(t,2), -mu*g*m)
eq

Let’s look at the properties of the differential equation:

sym.ode_order(eq, x(t))

sym.classify_ode(eq)

('factorable',
 'nth_algebraic',
 'nth_linear_constant_coeff_undetermined_coefficients',
 'nth_linear_euler_eq_nonhomogeneous_undetermined_coefficients',
 'nth_linear_constant_coeff_variation_of_parameters',
 'nth_linear_euler_eq_nonhomogeneous_variation_of_parameters',
 'nth_algebraic_Integral',
 'nth_linear_constant_coeff_variation_of_parameters_Integral',
 'nth_linear_euler_eq_nonhomogeneous_variation_of_parameters_Integral',
 '2nd_nonlinear_autonomous_conserved',
 '2nd_nonlinear_autonomous_conserved_Integral')

Let’s solve it:

solution = sym.dsolve(eq, x(t),
                     ics={
                         x(0):0,                     #at time 0s the displacement is zero
                         x(t).diff(t).subs(t, 0): v0 #at time 0s the velocity equals v0
                     })
solution

We retrieve the right-hand side of the equation like this (rhs - right-hand side):

solution.rhs

Additional

sympy.mechanics

sympy has built-in support for classical mechanics (documentation). A comprehensive tutorial was presented at the scientific conference SciPy 2016.

---

Review questions

---

Basics of object-oriented programming

Question 1: Prepare a class ComplexNumber whose constructor stores two given arguments in the attributes self.re and self.im. Create an object of the prepared class and print re and im.

Question 2: Extend the ComplexNumber class with a method conjugate() that enables complex conjugation of objects.

Question 3: Extend the class from the previous problem with a method absolute() that computes and returns the absolute value of the complex number, \(|a + \text{i}\,b| = \sqrt{a^2 + b^2}\).

Question 4: Prepare a class RealNumber that inherits from the ComplexNumber class. In the __init__ function it should call the constructor of the main class, assign the given value to the real part, and always assign the value 0 to the imaginary part.

---

SymPy

Question 5: Using the SymPy package, define the necessary symbols and define the expressions for the kinetic and potential energy \(E_k\), \(E_p\), and the Lagrangian \(L\) of a mass-spring system. Define all symbols as positive real numbers.

\[E_k = \frac{m~\dot{x}(t)^2}{2}\]

\[E_p = \frac{k~x(t)^2}{2}\]

\[L = E_k - E_p\]

Question 6: Using the equation written below, set up the equation of motion of the mass-spring system.

\[\frac{d}{dt}\Big(\frac{\partial L}{\partial \dot{x}(t)}\Big) - \frac{\partial L}{\partial x(t)} = 0\]

Question 7: Using a SymPy function, solve the obtained differential equation. Define a new symbol w and in the obtained solution perform the substitution: \(\sqrt{\frac{k}{m}} = \omega\). Store the right-hand side of the equation in the variable displacement (use solution.args).

Question 8: Determine the values of the integration constants 'C1', 'C2' by substituting the initial conditions for displacement and velocity into the obtained expression for the displacement:

\[ x~(t=0) = 1\]

\[ v~(t=0) = \dot{x}~(t=0) = 0\]

#Initial conditions:
x0 = 1
v0 = 0

Question 9: Substitute the obtained values of the constants 'C1' and 'C2' into the previously determined expression displacement.

Plot the obtained dependence \(x(t)\) over the interval \(0 \leq t \leq 2\) s, replacing the parameter w with the value \(2\pi\).

Question 10: Determine, to 5-place accuracy, the time t0 when the system under consideration is for the first time at the equilibrium position (\(x(t_0) = 0)\).

Question 11: Based on the final expression for the displacement \(x(t)\), create a numpy function to compute the displacement of the mass from the equilibrium position as a function of time (\(\omega = 2\pi\)).

Plot the result over the interval \(0 \leq t \leq 2\) s using the matplotlib library.