This is Words and Buttons Online — a collection of interactive #tutorials, #demos, and #quizzes about #mathematics, #algorithms and #programming.

SymPy makes math fun again

I remember my own struggle with calculus at university. Limits, integrals, differential equations. Lots of practice, lots of homework. Pages and pages of exercises. I loved math, loved the connection between algebra and geometry, loved the very pleasure of solving problems by making different concepts work together. But I hated doing the “paperwork”.

Taking it seriously, I still studied through the semester, studied harder the week before the exam, studied even harder the night before. I got 62/100. That's 1 point above the lowest possible passing grade.

Well, maybe math is not for everyone. But wait a minute! The next semester I took part in the Math Olympiad, went through the faculty round, then through the university round, went to the nationals, and even managed to score a few points there. Which counted as a pass on that semester's exam.

My professors were proud of me. And for almost a year, they thought that the first semester's score was a mistake. Until in the third semester, I scored 65/100.

¯\_(ツ)_/¯

Mathematics is a lot of things. It's the fun of problem-solving, it's the excitement of discoveries, it's the pride of accomplishments, and it's a ton of tedious computations, too. I never liked that last part. I still don't. That's why I'm so happy to live in the XXI century since I can give it away to computers and still enjoy the first three.

SymPy it is

According to the official site:

SymPy is 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.

Symbolic mathematics? Algebra system? Sounds exclusive. Sounds like it's a tool for practicing mathematicians and maybe the students who desperately want to become one. But it isn't. It is for practicing engineers who have just enough math knowledge to state a problem but not enough patience to solve it.

It does integration and differentiation. It finds limits, and it expands power series. It solves equations and systems of equations. It can even do statistics. And it does it all just like you would do on a math exam yourself. It doesn't just compute numbers, it computes formulas.

It is a Python library that does the boring part of math for you. Moreover, it does the math fast, accurately, and without angst. In other words, it is everything I'm not.

Let SymPy do some math for you

Let's say we want to model the sine function with a polynomial. Let's pretend we have a reason. Now how do we do that?

How about, we gather all we know about the sine function in one bowl and let SymPy do the rest? Sounds good?

sin(x)

We know that sin(0) = 0, right? Everybody knows that. Also, as you might have heard, you can approximate a small sine with its own argument: sin(0.001) ≃ 0.001. This means that the derivative of sine in 0 is 1.

The sine function climbs from 0 to π/2 and then it starts going down. In π/2 it's neither climbing nor descending, it's at its peak. This means that the derivative of sin(π/2)' = 0. Also, since the sine only climbs to 1, sin(π/2) = 1.

The sine is symmetrical relative to π/2. This means that it descends exactly as it climbs, and at the point where it reaches full π, its derivative is sin(π)' = -1, and the function itself sin(π) = 0.

Also, one less known (but easily computable) fact, the integral of the sine from 0 to π/2 is 1. Let's throw this into the bowl as well.

So, we have 7 facts. This implies 7 equations. This implies that our polynomial will have 7 coefficients and the highest degree will be then 6.

Let's translate it all into Python.

from sympy import *

a, b, c, d, e, f, g, x = symbols('a b c d e f g x')

sine = a*x**6 + b*x**5 + c*x**4 + d*x**3 + e*x**2 + f*x + g
sine_d = diff(sine, x)
sine_i = integrate(sine, x)

the_system = [
    sine_i.subs(x, pi / 2) - sine_i.subs(x, 0) - 1,
    sine_d.subs(x, 0) - 1,
    sine_d.subs(x, pi / 2),
    sine_d.subs(x, pi) + 1,
    sine.subs(x, 0),
    sine.subs(x, pi / 2) - 1,
    sine.subs(x, pi)
]

res = solve(the_system, (a, b, c, d, e, f, g))

print(res)

You can browse the tutorial if you like but it should be more or less clear as it is. sine is our polynomial model. diff is differentiation. integrate is integration. solve is solve. Simple!

Our solution now is a Python dictionary of formulas:

{a: (-448*pi - 28*pi**2 + 1680)/pi**7,
 b: (-5040 + 84*pi**2 + 1344*pi)/pi**6,
 c: (-1440*pi - 95*pi**2 + 5460)/pi**5,
 d: (-2520 + 50*pi**2 + 640*pi)/pi**4,
 e: (-96*pi - 12*pi**2 + 420)/pi**3,
 f: 1,
 g: 0}

Which is nice but not be exactly what we expected. We need the coefficients for a polynomial model and SymPy gave us the way to compute these coefficients from π instead. Well, it's what it does, it computes things symbolically not numerically. Or does it?

Look! I'll add one more line and SymPy will compute our coefficients as floating-point numbers.

from sympy import *
from math import pi

a, b, c, d, e, f, g, x = symbols('a b c d e f g x')

sine = a*x**6 + b*x**5 + c*x**4 + d*x**3 + e*x**2 + f*x + g
sine_d = diff(sine, x)
sine_i = integrate(sine, x)

the_system = [
    sine_i.subs(x, pi / 2) - sine_i.subs(x, 0) - 1,
    sine_d.subs(x, 0) - 1,
    sine_d.subs(x, pi / 2),
    sine_d.subs(x, pi) + 1,
    sine.subs(x, 0),
    sine.subs(x, pi / 2) - 1,
    sine.subs(x, pi)
]

res = solve(the_system, (a, b, c, d, e, f, g))

for var, exp in res.items():
    print(var, srepr(exp))

Have you noticed which line it is? Anyway, the result is now this:

{a: -0.00125233934372311,
 b: 0.0118030202461258,
 c: -0.00492072682786551,
 d: -0.163234062440025,
 e: -0.000907801926098641,
 f: 1.00000000000000,
 g: 0.0}

The line was from math import pi. This overloads the pi in the scope to be a floating-point number, not a symbol. And voilà — SymPy is now numeric. We can take these numbers, put them into our polynomial and it will be our model.

The polynomial model put over original sin(x)

The model works wonders in the [0; π] range. Outside this range, the model diverges from the sine but we never specified that it shouldn't.

The polynomial bears all the properties we told SymPy about with our equations and nothing more.

A shameless plug.

In my Geometry for Programers book, I used a less graphic but more practical example. I modeled the sine on [0, π/2] range with only four non-zero coefficients and then compared my model to a conventional one obtained as a power series. Mine appeared to be an order of magnitude more precise.

Now let SymPy write some code for us

Let's say we want something that looks like the sine but isn't. Something we can tweak interactively. Like the sine but with the movable middle point.

We can start the same way. Let's retain the point locations. The integral criterion and the derivatives in the endpoints criteria are better to be lifted. But that's just my opinion, you can try to retain them yourself and see what'll happen.

Anyway, the code that states the problem now looks like this:

from sympy import *

a, b, c, d, x, px, py = symbols('a b c d x px py')

sine = a*x**3 + b*x**2 + c*x + d
sine_d = diff(sine, x)

the_system = [
    sine_d.subs(x, px),
    sine.subs(x, 0),
    sine.subs(x, px) - py,
    sine.subs(x, pi)
]

res = solve(the_system, (a, b, c, d))

print(res)

When run, it prints out that:

{a: (-2*px*py + pi*py)/(px**4 - 2*pi*px**3 + pi**2*px**2),
 b: (3*px**2*py - pi**2*py)/(px**4 - 2*pi*px**3 + pi**2*px**2),
 c: (-3*pi*px*py + 2*pi**2*py)/(px**3 - 2*pi*px**2 + pi**2*px),
 d: 0}

This is already code-like, you can copy-paste this to your program and maul it a bit so it would fit the syntax of your favorite language. Or you can ask SymPy to do that for you as well!

Just replace print(res) with print(jscode(res)) and SymPy will write the solution in JavaScript!

{a: (-2*px*py + Math.PI*py)/(Math.pow(px, 4) -
    2*Math.PI*Math.pow(px, 3) +
    Math.pow(Math.PI, 2)*Math.pow(px, 2)),
 b: (3*Math.pow(px, 2)*py - Math.pow(Math.PI, 2)*py)/
    (Math.pow(px, 4) - 2*Math.PI*Math.pow(px, 3) +
    Math.pow(Math.PI, 2)*Math.pow(px, 2)),
 c: (-3*Math.PI*px*py + 2*Math.pow(Math.PI, 2)*py)/
    (Math.pow(px, 3) - 2*Math.PI*Math.pow(px, 2) +
    Math.pow(Math.PI, 2)*px),
 d: 0}

The formulas are now completely runnable code. Here I used them to make this very interactive widget ↓

You can drag the point. SymPy-generated code will recompute the curve.

SymPy can also write code in Rust, C++, Fortran, Matematica, and more. Most of the non-esoterical languages are already there.

They say that AI will replace programmers in the future. Well, for a certain part, the future is here. SymPy already writes your code for you and not just some boilerplate. It really solves a mathematical problem and turns the solution into something you can run. Well, if writing code means solving equations without any creative thought involved, I don't mind being replaced.

Conclusion

I disagree with the “math is not for everyone”. Math is vast and diverse. Surely there is enough math for everyone. But some parts of it are better left to computers.

Hope this demonstration shows that computers are happy to help you. It's not that hard to state a problem and it's even easier to interpret results. And it's fun.

If my humble introduction to SymPy will help you enjoy math long after your final calculus exam, it'll make me happy too.