Symbolic Algebra

Pystencils lifts numerical programming to the level of symbolic algebra. Instead of writing code and interacting with numerical data directly, we declaratively describe a numerical kernel by constructing and manipulating mathematical expressions. This allows us to, among other things,

• define numerical methods on an algebraic level of abstraction while hiding the intricate implementation details of high-performance computing;

• derive numerical rules from first mathematical principles semi-automatically, using tools like symbolic differentation or integration;

• and perform extensive algebraic simplifications informed by domain-specific knowledge.

To make all of this possible, pystencils relies on the versatile Python package SymPy. SymPy is a computer algebra system implemented as a Python library. Through it, we can construct and rewrite mathematical expressions as objects. As SymPy is the very foundation of pystencils, we will begin this lesson with a brief introduction into its most basic features. Throughout this tutorial suite, when using features of SymPy, we will introduce them shortly. If you are going to keep using pystencils productively, however, we strongly recommend you take the time to work through the official SymPy tutorial.

Getting Started with SymPy¶

To begin using SymPy, we first need to import it. SymPy usually goes by the canonical abbreviation sp:

import sympy as sp

Let’s start by considering the basic building blocks of algebraic expressions: numbers and symbols. Numbers in SymPy are distinct mathematical entities, and are always represented exactly. Creating a constant SymPy expression is very different from using a floating-point literal or computing a constant numerically; while the former will-- at most --reduce to a canonical form through algebraic rules, the latter will evaluate to a floating-point value, introducing discretization and round-off errors.

Numbers, Algebraically¶

Let’s look at this in practice. The following cell computes the value of the fraction $\frac{1}{12}$ numerically:

1.0 / 12.0

The next cell, on the other hand, defines it algebraically, using sp.Rational:

sp.Rational(1, 12)

As you can see, the numerical value is inexact - after all, $\frac{1}{12}$ cannot be represented with only finitely many binary digits. The symbolic version, being stored as a fraction of integers, will always be mathematically accurate. This becomes more apparent when we add two fractions together:

1.0 / 12.0 + 3.0 / 31.0

sp.Rational(1, 12) + sp.Rational(3, 31)

You should see the result as a fraction, reduced to the common denominator $12 \cdot 31 = 372$. This is one of the powers of symbolic computing: Expressions can be algebraically simplified to their most efficient form without inducing rounding errors.

Symbols and Complex Expressions¶

Let’s consider the second fundamental building block: Symbols. While appearing similar, symbols in symbolic algebra are very different from variables in classical programming. A variable is a name for a memory location, referring a piece of data stored there. Symbols, on the other hand. are named placeholders for mathematical objects from some underlying domain. In SymPy, this domain, by default, are the real numbers $\mathbb{R}$.

Symbols are objects; they can be created using the sp.symbols factory:

x = sp.symbols("x")

You can create multiple symbols in one line, as well:

y, z, w  = sp.symbols("y, z, w")

Now, let’s assemble symbols and constants into more complex expressions. SymPy understands all the fundamental arithmetic operators of Python; when applied to expressions, they result in the construction of larger expressions from the operands. Thus, expressions, too, are objects. Here is an example:

expr = x + y * w - z**3
expr

Expressions are a tree structure: their inner nodes are operators and functions, while symbols and constants come in at the leaves. To see this in a bit more detail, we can use graphviz to visualize the expression tree:

import graphviz
tree = graphviz.Source(sp.printing.dotprint(expr))
tree

The tree, as you should now be able to see, does not exactly correspond to the printed form of the expressions. For instance, subtractions are actually internally modelled as a combination of addition and multiplication with -1.

# 1)
5 * x**2 - 3 * x + 2

# 2)
a, b, c = sp.symbols("a, b, c")
a * x**2 + b * x + c

# 3)
sp.exp(x) * (sp.sin(x) + sp.cos(x))

def poly1d(coeffs):
    x = sp.symbols("x")
    return sum(c * x**i for i, c in enumerate(coeffs))

So much for the basics of symbolic algebra. There is, of course, much more to explore here. SymPy allows us to automatically simplify and rewrite expressions according to the rules of algebra. We can compute derivatives and integrals, perform linear algebra with vectors and matrices, and much more. To dive deeper into these things, we refer to the SymPy tutorial.

This concludes the lesson on symbolic algebra; You may now continue with the next lesson, where we introduce the first, and most important, concept which pystencils adds to SymPy: fields.