3.4 Basics of python

“By relieving the brain of all unnecessary work, a good notation sets it free to concentrate on more advanced problems.”

— Alfred North Whitehead, An Introduction to Mathematics, 1911

This chapter covers the Python tools we use throughout this book: SymPy for symbolic mathematics, NumPy for numerical computation, and Matplotlib for visualization. The focus here is on the subset needed for mechanics. For a comprehensive treatment of computational thinking, applied mathematics and Python applications, see python.ju.se.

We also use MechanicsKit, a lightweight package that provides convenience functions for engineering education: LaTeX rendering of arrays, MATLAB-style plotting of symbolic expressions with fplot, 1-based indexing for FEM workflows, and more. Install it with:

uv pip install git+https://github.com/cenmir/MechanicsKit.git

Working in a notebook

The code in this book is written to be run in a notebook. A notebook is a file made of cells. A text cell holds prose and equations, and a code cell holds a few lines of Python. Running a code cell executes its lines in order and shows the result directly below it, so a calculation and its outcome stay together on the page. The names a cell creates survive after it has run, and every later cell can use them, so the chapters build a problem up cell by cell instead of in one long program. Two rules of the notebook matter from the first cell on. Only the value of the last line of a cell is displayed, and nothing from the lines above it. And a cell sees whatever has been run before it, not whatever stands above it, so running cells out of order can leave a name holding an old value; when in doubt, run the notebook from the top.

There are two notebook formats to choose from and this book supports both. A Jupyter notebook, a file ending in .ipynb, is edited in VS Code or in the browser at jupyter.ju.se (for Jönköping University students), and it is the format these chapters are written in. A marimo notebook is a plain .py file that opens in the browser. It tracks which cells depend on which names and reruns them for us when a value changes, so the out-of-order problem disappears, at the price that a name may be defined in only one cell. Both run the same Python and the same packages, and every code cell in this book can be pasted into either. The installation chapter sets both up.

Symbolics using SymPy

To get started with symbolic computation in Python, we will use the SymPy library. SymPy is a Python library for symbolic mathematics, also known as computer algebra system (CAS). We can use SymPy to perform algebraic manipulations, calculus, solve equations, and much more. We load the SymPy library using the following command. import loads the library, and as sp gives it a short name for the rest of the notebook, so that every SymPy function is reached as sp. followed by its name:

import sympy as sp

Let us showcase some of the basic functionalities of SymPy using an example, we will solve the equation \(x^2 = 4\). SymPy has to be told that \(x\) stands for an unknown, since until now every x in a program would have been a number, and that is what the first line arranges.

The first line, x = sp.symbols('x'), is the first Python statement in the book, so we read it in parts. The name x on the left is ours to choose. The = is not an equality: it stores whatever the right side produces under that name. On the right, sp is the short name we gave SymPy in the import, the dot reaches into it, and symbols is a function found there. The parentheses call the function, and 'x' in quotes is a string, the text SymPy will print for the symbol. The two x are unrelated, one a Python name and the other a label, and they could differ. Since = is taken, an equation needs its own function, Eq, with the two sides as its arguments, and ** raises to a power. Finally solve takes the equation and the unknown and returns every value of the unknown that satisfies it.

1x = sp.symbols('x')
2eq = sp.Eq(x**2, 4)
3sol = sp.solve(eq, x)
4sol

1: Make an unknown named x.
2: Build the equation as an object. Python’s = assigns, so an equation needs Eq.
3: Ask for every value of x that satisfies eq. The answer is a list.
4: A name alone on the last line of a cell is displayed. The lines above it display nothing.

[-2, 2]

Three names now hold three kinds of object. Before running the next cell, guess what each of them holds. The function type reports it:

type(x), type(eq), type(sol)

(sympy.core.symbol.Symbol, sympy.core.relational.Equality, list)

So x is a SymPy Symbol, eq an Equality and sol a plain Python list. Displayed, the list is readable for two numbers, but not once the results are vectors and long expressions, so MechanicsKit provides two display tools. The first is the pipe | la, which renders whatever stands to its left as LaTeX:

1from mechanicskit import la
2sol | la

1: This form imports one name out of a package, so we write la and not mechanicskit.la.
2: Between numbers | means “or”. MechanicsKit reuses it as a pipe: the object on the left is handed to la on the right, and the display is the value of this last line.

\[ {\def\arraystretch{1.0}\begin{bmatrix} -2 \\ 2 \end{bmatrix}} \]

la needs no further input, and that makes it the quick way to look at a result while working a problem out. What it shows is anonymous: a column of numbers with no name attached. The second tool, ltx, lets us write the display ourselves. Its arguments alternate between LaTeX strings and objects, and it joins them into one line: here the string x_1 =, the first root, the string ,\quad x_2 = and the second root. The strings are raw strings, r"...", so that backslashes such as the one in \quad, a LaTeX space, reach LaTeX unchanged:

from mechanicskit import ltx
1ltx(r"x_1 = ", sol[0],
2    r",\quad x_2 = ", sol[1])

1: A LaTeX string and then the object it labels. sol[0] is the first root, since Python counts from zero.
2: The next string starts with a comma and \quad, a LaTeX space, so the second root follows on the same line. sol[1] is the second root.

\[ x_1 = -2,\quad x_2 = 2 \]

We can also attach assumptions to a symbol. If we know that \(x\) is real, we create it with real=True, and if we know that it is positive, with positive=True. This is a keyword argument: the name says which setting is meant and True switches it on, and it follows the ordinary arguments in the call. solve then discards every root that violates the assumption, and one string and one object are enough for ltx:

1x = sp.symbols('x', positive=True)
eq = sp.Eq(x**2, 4)
2sol = sp.solve(eq, x)
3ltx(r"x = ", sol[0])

1: The name x is reused for a new symbol that carries the assumption; the old symbol is forgotten with it, and so is the old eq, so it is built again.
2: solve honours the assumption: \(-2\) is not positive, so the list holds one root.
3: With one root, item 0 is the whole answer.

\[ x = 2 \]

Example 1 - Creating a vector and equation system

We create a force vector \(\bm F\) whose components \(F_x\) and \(F_y\) are unknown, with all forces in newtons. One call of symbols makes as many symbols as the string names, separated by spaces, and the two names on the left receive them in order. The doubled letter in FF is our convention for a name that holds a vector, as \(\bm F\) is bold on paper. The objects passed to ltx can be matrices as well as numbers, and a SymPy Matrix is displayed as a column vector. We build \(\bm F\) as a sp.Matrix and not as a plain Python list; Lists and matrices at the end of the chapter explains why. In the string, \bm sets the name in bold, the book’s notation for a vector:

1F_x, F_y = sp.symbols('F_x F_y')
2FF = sp.Matrix([F_x, F_y, 0])
ltx(r"\bm F = ", FF)

1: Two unknowns from one call. The underscore in a label becomes a subscript when the symbol is printed.
2: The list in brackets supplies the components. The third is 0 because the book keeps every vector three-dimensional, so that the cross product is always available.

\[ \bm F = \left[\begin{matrix}F_{x}\\F_{y}\\0\end{matrix}\right] \]

A second force \(\bm F_1\) of 250 N acts at 38° to the \(x\) axis, and the two forces together must give a resultant of 1000 N along \(x\). Magnitude times direction gives \(\bm F_1\), and the resultant condition is a vector equation,

\[ \bm F_1 = 250\begin{bmatrix} \cos 38^\circ \\ \sin 38^\circ \\ 0 \end{bmatrix}, \qquad \bm F + \bm F_1 = \begin{bmatrix} 1000 \\ 0 \\ 0 \end{bmatrix} \]

SymPy works in radians, so we first define sind and cosd that take degrees. A function is defined with def, a name and its parameters in parentheses; the indented lines below form its body, and return hands the result back to the caller. The x inside sind is a placeholder for whatever is passed in and has nothing to do with the symbol x outside. The equation eq has a left side eq.lhs and a right side eq.rhs, and we pass both to ltx with the string = between them:

from sympy import rad

1def sind(x):
2    return sp.sin(rad(x))
def cosd(x):
    return sp.cos(rad(x))

3FF1 = 250*sp.Matrix([cosd(38), sind(38), 0])

4eq = sp.Eq(FF + FF1, sp.Matrix([1000, 0, 0]))
5ltx(r"\bm F + \bm F_1 = ", eq.lhs,
6    r" = ", eq.rhs)

1: The problem states degrees and SymPy’s sin wants radians, so the conversion is written once, here, and not at every use.
2: rad converts, sp.sin acts on the result, and what follows return is what a call of sind gives back.
3: Magnitude times unit direction, \(F_1 \bm e_1\). The scalar multiplies every component.
4: Eq on two matrices is one equation per component, three at once.
5: eq keeps both sides, and .lhs fetches the left one: the sum, with the unknowns still in it.
6: .rhs fetches the right one, the required resultant. Displaying both shows the equation as SymPy holds it.

\[ \bm F + \bm F_1 = \left[\begin{matrix}F_{x} + 250 \cos{\left(\dfrac{19 \pi}{90} \right)}\\F_{y} + 250 \sin{\left(\dfrac{19 \pi}{90} \right)}\\0\end{matrix}\right] = \left[\begin{matrix}1000\\0\\0\end{matrix}\right] \]

Two unknowns need two equations, and the vector equation supplies three, of which the third, \(0 = 0\), holds by itself. Solving gives the two unknown components.

1sol = sp.solve(eq, [F_x, F_y])

1: The list names the unknowns to solve for. With one unknown the symbol alone was enough; with two or more they go in brackets, and solve uses the component equations together.

The result is a dictionary that maps each unknown to its solution. Piped through la, a dictionary becomes an aligned system of equations, one line per unknown:

sol | la

\[ \begin{aligned} F_{x} &= 1000 - 250 \cos{\left(\dfrac{19 \pi}{90} \right)} \\[9.0pt] F_{y} &= - 250 \sin{\left(\dfrac{19 \pi}{90} \right)} \end{aligned} \]

Each solution is reached with its symbol as the key, as in sol[F_x]. The solutions are exact SymPy expressions that still contain \(\pi\). They can enter further calculations, or be evaluated to floating point numbers with evalf. For the numeric result, the precision argument of ltx sets the number of significant figures and a final string adds the unit:

1FF_num = sp.Matrix([sol[F_x], sol[F_y], 0]).evalf()
2ltx(r"\bm F = ", FF_num, r"\ \text{N}", precision=4)

1: The dictionary is read with the symbols as keys, the two values are assembled back into a vector, and evalf turns every entry numeric in one go. A dot after a closing parenthesis applies the method to what the call just built.
2: precision rounds the display only; FF_num keeps its full precision for any later calculation.

\[ \bm F = \left[\begin{matrix}803.0\\-153.9\\0\end{matrix}\right]\ \text{N} \]

A report often wants the exact result and its value side by side, and that layout is ours to choose with ltx. Passing aligned=True lets a display span several lines: \\ starts a new line, and & marks the point at which the lines are aligned, here the equals sign. evalf(4) rounds each value to four significant figures. The call itself spans two lines, which Python allows inside an open parenthesis, and the two strings that meet at the line break are joined into one:

1ltx(r"F_x &= ", sol[F_x], r"\approx ", sol[F_x].evalf(4), r"\ \text{N} \\ "
2    r"F_y &= ", sol[F_y], r"\approx ", sol[F_y].evalf(4), r"\ \text{N}",
3    aligned=True)

1: One display line: the exact value, then its four-figure evaluation, then the unit. The \\ closing the string ends the line, and the string continues on the next code line.
2: The second display line, built the same way. Both hold an & just before the equals sign.
3: aligned=True turns the & marks into alignment points; without it they would be an error.

\[ \begin{aligned}F_x &= 1000 - 250 \cos{\left(\dfrac{19 \pi}{90} \right)}\approx 803.0\ \text{N} \\ F_y &= - 250 \sin{\left(\dfrac{19 \pi}{90} \right)}\approx -153.9\ \text{N}\end{aligned} \]

Lists and matrices

A Python list is a general container. It holds items of any kind in order, and its operators act on the container, not on the items: adding two lists joins them, and multiplying a list by an integer repeats it.

a = [1, 2, 0]
b = [3, 4, 0]
1a + b, 2*a

1: Two results on the last line, separated by a comma, are displayed together as a pair.

([1, 2, 0, 3, 4, 0], [1, 2, 0, 1, 2, 0])

Neither result is what a vector should give. A sp.Matrix follows the rules of linear algebra, so addition works component by component and a scalar multiplies every component:

aa = sp.Matrix([1, 2, 0])
bb = sp.Matrix([3, 4, 0])
1ltx(r"\bm a + \bm b = ", aa + bb,
2    r",\quad 2\bm a = ", 2*aa)

1: + on two matrices adds component by component, the vector sum.
2: A number times a matrix scales every component. The string opens with a comma and \quad so the two results share one line, separated by a space.

\[ \bm a + \bm b = \left[\begin{matrix}4\\6\\0\end{matrix}\right],\quad 2\bm a = \left[\begin{matrix}2\\4\\0\end{matrix}\right] \]

The vector operations of mechanics are methods of the matrix: dot for the scalar product, cross for the vector product and norm for the length,

\[ \bm a \cdot \bm b = a_x b_x + a_y b_y + a_z b_z, \qquad \bm a \times \bm b = \begin{bmatrix} a_y b_z - a_z b_y \\ a_z b_x - a_x b_z \\ a_x b_y - a_y b_x \end{bmatrix}, \qquad |\bm a| = \sqrt{\bm a \cdot \bm a} \]

A method is a function that belongs to the object. It is called with a dot after the object’s name, and the other vector, when one is needed, goes in the parentheses. The three results go into one ltx call that spans three lines, one per operation:

1ltx(r"\bm a \cdot \bm b = ", aa.dot(bb),
2    r",\quad \bm a \times \bm b = ", aa.cross(bb),
3    r",\quad |\bm a| = ", aa.norm())

1: The scalar product of aa with bb is a number, so it is displayed inline.
2: The vector product is a vector, so it is displayed as a column. The order matters: bb.cross(aa) would give the opposite sign.
3: The length involves aa alone, so norm takes no argument. The empty parentheses still make the call.

\[ \bm a \cdot \bm b = 11,\quad \bm a \times \bm b = \left[\begin{matrix}0\\0\\-2\end{matrix}\right],\quad |\bm a| = \sqrt{5} \]

In sp.Matrix([1, 2, 0]) the list only supplies the components, and the matrix it builds is a column vector. A matrix can also hold symbols, and that is what Example 1 relies on: \(\bm F\) has unknown components, adding \(\bm F_1\) to it gives a vector equation, and solve finds the unknowns from all of its components at once. NumPy arrays share the arithmetic but hold numbers only, and we use them once a problem has become numerical.