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, the connection between algebra and geometry, 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 throught the university round, went to the national round and scored quite high there. Which counted as a pass on that semester's exam.

Professors were proud of me. They thought that the first semester's score was a mistake. The third semester I scored 65/100.

¯\_(ツ)_/¯

Mathematics is a lot of things. It's fun of problem-solving, it's an excitement of discoveries, it's a pride of accomplishments, and it's a ton of tedious computations, too. I never liked the last part. Couldn't make it right. 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 something for practicing mathematicians only. But it isn't. It's 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 series. It solves equations and systems of equations. It can even do statistics. And it does it all just like you would do on an exam. It doesn't just compute numbers, it computes formulas.

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

Let SymPy do some math for you

Let's say you want to model the sine function with a polynomial. Let's pretend you 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 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 or descending, it's in 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 system 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))

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

You can browse the tutorial if you like but it seems more or less clear as it is. sine is our polynomial model. diff is differentiation. integrate is integration. solve is solve. Simple! The only cryptic thing here is srepr. It turns the solution into a special form of expression that looks like this:

  (f, 'Integer(1)')
  (c, 'Mul(Integer(5), Pow(pi, Integer(-5)), Add(Mul(Integer(-1), Integer(288), pi), Mul(Integer(-1), Integer(19), Pow(pi, Integer(2))), Integer(1092)))')
  (d, 'Mul(Integer(10), Pow(pi, Integer(-4)), Add(Integer(-252), Mul(Integer(5), Pow(pi, Integer(2))), Mul(Integer(64), pi)))')
  (a, 'Mul(Integer(28), Pow(pi, Integer(-7)), Add(Mul(Integer(-1), Integer(16), pi), Mul(Integer(-1), Pow(pi, Integer(2))), Integer(60)))')
  (e, 'Mul(Integer(12), Pow(pi, Integer(-3)), Add(Mul(Integer(-1), Integer(8), pi), Mul(Integer(-1), Pow(pi, Integer(2))), Integer(35)))')
  (g, 'Integer(0)')
  (b, 'Mul(Integer(84), Pow(pi, Integer(-6)), Add(Integer(-60), Pow(pi, Integer(2)), Mul(Integer(16), pi)))')

Ok, this is not really helpful. We need the coefficients and SymPy gave us the way to compute these coefficients from π. 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:

(f, "Float('1.0', precision=53)")
(c, "Float('-0.0049207268278655082', precision=53)")
(d, "Float('-0.16323406244002545', precision=53)")
(a, "Float('-0.0012523393437231091', precision=53)")
(e, "Float('-0.00090780192609864056', precision=53)")
(g, "Float('0.0', precision=53)")
(b, "Float('0.011803020246125803', precision=53)")

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.

Polynomial model put over original sin(x)

The model works wonders in the [0; π] range. Outside this range it 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.

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 along with the derivatives in the endpoints are better to be lifted. But that's just my opinion, you can try to retain them yourself and see what'll happen.

The code that states the problem look 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))

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

When run, it prints out that:

  (c, "Mul(pi, Pow(Symbol('px'), Integer(-1)), Symbol('py'), Add(Mul(Integer(-1), Integer(3), Symbol('px')), Mul(Integer(2), pi)), Pow(Add(Pow(Symbol('px'), Integer(2)), Mul(Integer(-1), Integer(2), pi, Symbol('px')), Pow(pi, Integer(2))), Integer(-1)))")
  (b, "Mul(Pow(Symbol('px'), Integer(-2)), Symbol('py'), Add(Mul(Integer(3), Pow(Symbol('px'), Integer(2))), Mul(Integer(-1), Pow(pi, Integer(2)))), Pow(Add(Pow(Symbol('px'), Integer(2)), Mul(Integer(-1), Integer(2), pi, Symbol('px')), Pow(pi, Integer(2))), Integer(-1)))")
  (a, "Mul(Pow(Symbol('px'), Integer(-2)), Symbol('py'), Add(pi, Mul(Integer(-1), Integer(2), Symbol('px'))), Pow(Add(Pow(Symbol('px'), Integer(2)), Mul(Integer(-1), Integer(2), pi, Symbol('px')), Pow(pi, Integer(2))), Integer(-1)))")
  (d, 'Integer(0)')

And now it's time to get back to that cryptic srepr thing from before. It stands for symbolic expression, and it's basically a ready-made syntax tree. It's like in Lisp only with somewhat more conventional parentheses. And since it's written in conventional C-style function form, it's incredibly easy to write an interpreter for it (see: the code on GitHub).

Here is a simple tool that turns SymPy s-expressions into runable JavaScript code.


Just press the button — and that's it! Of course, you can prettify the code however you like. But you can also use in raw like this (see: the example from this very page).

You can drag the point

But the best part of it, unless you're writing in some esoteric languages, you don't even need an interpreter. SymPy can print out JavaScript code from the box with a single jscode function! It can also print in Rust, C++, Fortran, Matematica, and more. It's already there.

Conclusion

I disagree with “math is not for everyone”. Math is vast and diverse. Surely there is enough math for everyone. But some parts of it are best 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.