Playing with SymPy « Reasonable Deviations

Playing with SymPy

By Rod Carvalho

I have been playing with SymPy lately. SymPy is, basically, a Python library for Symbolic Mathematics. It’s precisely the kind of CAS I was looking for: simple, lightweight, easy to use, syntactically tasteful, open-source, and free. SymPy just makes sense. Simplicity über alles.

Moreover, SymPy is based on Python, and I do love Python. It’s great to be able to use a CAS without having to learn yet another scripting language with a bastard syntax. There are dozens of CAS‘es out there, each with its own language. Why shouldn’t one be able to use a CAS using a “traditional” language such as C++ and Python? There’s SAGE, which also uses Python, but SAGE is anything but lightweight. There were a couple of CAS‘es based on C++ too, but they were kind of heavy too. Why not keep it simple?

SymPy was started in 2005, by Ondřej Čertík, a Physics student from Czech Republic. Kudos to Ondřej for his initiative! And thanks to all of those who contributed to the SymPy project! Here are some slides on SymPy, by Ondřej:

SymPy – Python library for symbolic mathematics (PDF – 223 KB)

You can download SymPy from here. If you don’t have Python installed in your computer, you can play with the SymPy online shell (which runs on the Google App Engine).

__________

An example

Let’s consider the system of polynomial equations in $\mathbb{R}[x_1, x_2, x_3]$

$\left\{\begin{array}{rl}x_1 + x_2 + x_3 = y_1 + y_2 + y_3 \\x_1^2 + x_2^2 + x_3^2 = y_1^2 + y_2^2 + y_3^2\\x_1^3 + x_2^3 + x_3^3 = y_1^3 + y_2^3 + y_3^3\end{array}\right.$ ,

where $y = (y_1, y_2, y_3) \in \mathbb{R}^3$ is a known vector (whose entries are distinct). We will denote $x = (x_1, x_2, x_3)$ . Let function $g_k: \mathbb{R}^3 \rightarrow \mathbb{R}$ be defined as

$g_k(z) = \displaystyle\sum_{i=1}^3 z_i^k$ .

The system of polynomial equations can thus be rewritten as

$\left\{\begin{array}{rl}p_1(x) = 0\\p_2(x) = 0\\p_3(x) = 0\end{array}\right.$

where $p_k(x) = g_k(x) - g_k(y)$ for all $k \in \{1,2,3\}$ . According to my conjecture, this system of equations should have exactly $3! = 6$ solutions. Solving the system is not very easy because the equations $p_k(x) = 0$ are nonlinear in variables $\{x_1, x_2, x_3\}$ . A CAS would be useful.

For example, we could implement function $g_k: \mathbb{R}^3 \rightarrow \mathbb{R}$ with the following Python script:


def g(z,k):
    """Computes g_k(z), where z is a list of reals and k is a positive integer"""

    # checks function arguments for errors
    if len(z)==0:
        return "ERROR: First argument must be a non-empty list of symbols!"
    if (type(k) != int) or (k < 1):
        return "ERROR: Second argument must be a positive integer!"

    # computes g_k(z) and returns it
    acc = 0
    for i in range(0,len(z)):
        acc += (z[i])**k
    return acc

We could then compute the Gröbner basis of the set of polynomials $\{p_1(x), p_2(x), p_3(x)\}$ , and the solutions of the system of polynomial equations for a given $y$ , say, $y = (1,2,3)$ :


from sympy import *

# number of equations and symbolic variables
n = 3               

# defines set S = {1, 2,..., n}
S = range(1,n+1)

# declares symbolic variables
x = [Symbol('x%d' % i) for i in S]

# initializes y vector
#y = [Symbol('y%d' % i) for i in S]
y = S

# defines system of polynomial equations in variables [x1, x2, x3]
P = [g(x,k) - g(y,k) for k in S]

# computes Groebner basis and solutions of the system of polynomials
GB = groebner(P, order='lex')
Sols = solve_system(P)

# prints results
print "Groebner basis: %s" % GB
print "Solutions: %s" % Sols
print "There are %d solutions." % len(Sols)

The output is:

Gröbner basis: the Gröbner basis is a set of three polynomials $\{b_1(x), b_2(x), b_3(x)\}$ , where

$\begin{array}{l} b_1(x) = -6 + x_1 + x_2 + x_3\\ b_2(x) = 11 - 6 x_2 - 6 x_3 + x_2 x_3 + x_2^2 + x_3^2\\ b_3(x) = -6 + 11 x_3 - 6 x_3^2 + x_3^3\end{array}$ .

Solutions: there are $3! = 6$ solutions

$\{(1, 2, 3), (1, 3, 2), (2, 1, 3), (2, 3, 1), (3, 1, 2), (3, 2, 1)\}$ .

As expected, the solutions of the system of polynomial equations are the $3! = 6$ permutations of the elements of $y = (1,2,3)$ . Of course, this does not prove the conjecture. All it proves is that my conjecture works for the particular case where $n=3$ and $y = (1,2,3)$ .

__________

2008 Google Summer of Code and SymPy

2007 Google Summer of Code and SymPy

SymPy/sympycore (pure Python) up to 5x faster than Maxima

You may take away my Mathematica license but you may never take my freedom

MATLAB to drop Maple for Symbolic Toolbox

Computer Algebra Systems in odd places #1

Tags: Computer Algebra Systems, Ondřej Čertík, Python, Symbolic Mathematics, SymPy

This entry was posted on July 5, 2008 at 12:58 am and is filed under Software. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.

6 Responses to “Playing with SymPy”

Alex Says:
June 13, 2008 at 9:08 am | Reply
I’ve used SAGE for a total of 5 minutes, but I do have a strong opinion on it. It’s difficult for me to believe that you can sucessfully paste the illusion of a uniform interface over several CASes disparate in goals, philosophies, and abilities. Maybe SAGE has done that, but I also have the more personal issue that if I were to use it, I’d like to be proficient in using the various backends it calls. It just wouldn’t feel right otherwise :)

You could try JACAL or MockMMA for short nontrivial (but probably confusingly LISPy) CAS codes; I think you’d like the feel of SymPy, and I bet the code’s readable; Maxima for a more industrial sized codebase; or wait forever and maybe the Axiom people will get around to documenting their code (I personally can’t wait for that).

Either way, you’ve got to learn a lot of hardcore math (boring boring algebra) to understand their algorithms.
Alex Says:
July 2, 2008 at 10:50 pm | Reply
Maple looks ugly, Mathematica is syntactically annoying … and those are the only two worth considering for general use. I settled on Mathematica– next fall I’ll actually be co-teaching a course on it, so over the summer I must become some kind of expert :)

My hope is reserved for Axiom. The core is written in LISP with most algorithms implemented in the system’s extension language. I’m not sure about the syntax of the system; I haven’t seriously used it ever, and my dabblings were long ago, but the idea of a category theoretical CAS seems very powerful. On top of that, I messed around with Aldor (the latest of its extension languages), and I found it fun. Hopefully they’ll get around to documenting the code so people can go in and rip off the algorithms one by one and incorporate them into their own CASes.

That’s probably the ideal situation: having a compendium of code that groups can use to roll their own CASes that fit their own design specifications.
Ondrej Says:
July 5, 2008 at 11:53 pm | Reply
Nice to hear that SymPy is useful. I’ve added your post to:

http://docs.sympy.org/outreach.html#blogs-news-magazines

and also cited it here:

http://groups.google.com/group/sympy/browse_thread/thread/e33d7d03d24f2a9c

If you have any suggestions for improvements, please let us know on our mailinglist. We are also looking for new users and developers if anyone is interested. :)
rod. Says:
July 6, 2008 at 4:28 am | Reply
If you know Python, you’ll probably think that in the piece of code that implements function $g_k(\cdot)$ , the following lines
```
# computes g_k(z) and returns it
acc = 0
for i in range(0,len(z)):
    acc += (z[i])**k
return acc
```
don’t look very “Pythonic” (at all). I know, I know. These lines of code look as though they were translated from C or C++ to Python (I even use the “+=” operator). This is kind of “plain vanilla”.

It would be much cooler to replace the lines above with a single line of code that does it all:
```
# computes g_k(z) and returns it
return sum([elem**k for elem in z])
```
This would work if the function’s first argument were a list of reals. However, if it’s a list of symbols, then function sum() does not work. I suppose the problem is that sum() uses an operator which does not work on symbols. Operator overloading would be necessary.
Shubhendu Trivedi Says:
July 6, 2008 at 1:57 pm | Reply
Yes, i have been learning Python as I told you once. And indeed Python makes life easier. Actually when you i said i am walking on air i was reminded of this. Here the chap is walking on air (metaphorically). I got interested in python largely because of SimPy (as a substitute for StarLogo) but after starting with it i started loving it!
jkwiens.com: Learning Python Says:
March 20, 2009 at 11:35 am | Reply
[...] playing around with SymMath which is a computer algebra system for Python. If it wasn’t for Rod, I never would have known [...]

Reasonable Deviations

Playing with SymPy

6 Responses to “Playing with SymPy”

Leave a Reply