[Tutor] Python and Symbolic Math for beginners
Oscar Benjamin
oscar.j.benjamin at gmail.com
Mon Jun 17 17:13:49 CEST 2013
On 15 June 2013 10:53, Amit Saha <amitsaha.in at gmail.com> wrote:
> Would any of you have any teaching (or substantial self learning)
> experience with a library for Symbolic math?
I have taught using Maple. I haven't with Sympy but I do use it myself.
> I am currently exploring sympy (http://sympy.org) as part of writing a
> book chapter and would like to know if there any better/easier option
> out there which can successfully introduce symbolic math to young
> programmers.
In terms of free software sympy/sage are probably the best options. I
guess that my preference for teaching symbolic computation would be to
begin with isympy (sympy mode in the IPython console). I think that
this is initially better than the notebook modes of sympy, Maple, etc.
because it is clear what is input and what is output and what
functions are being called. That way you can gain a sense of the nuts
and bolts for how symbolic computation works without getting confused
by the "magic" of most of these systems. Later most people would
probably find it more practical to use the isympy notebooks or sage
for serious work.
If you use isympy with the -a flag then missing symbols are
automatically created solving the inconvenience that Jim referred to.
Also if you have a decent console you can pass '-p unicode' to get all
the Greek letters (this is automatic on e.g. Ubuntu but not on
Windows):
$ isympy -p unicode -a
IPython console for SymPy 0.7.2 (Python 2.7.5-32-bit) (ground types: python)
These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
Documentation can be found at http://www.sympy.org
In [1]: alpha + kappa
Out[1]: α + κ
Unfortunately the auto-symbols feature doesn't work so well when a
particular symbol is unexpectedly defined e.g. the beta function:
In [2]: alpha + beta
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-2-92b8e9c3b3ff> in <module>()
----> 1 alpha + beta
TypeError: unsupported operand type(s) for +: 'Symbol' and 'function'
Otherwise I think that sympy has a good feature set. Some examples:
Get the general formula for the sum of i**2 with i from 1 to n:
In [1]: summation(i**2, (i, 1, n))
Out[1]:
3 2
n n n
-- + -- + -
3 2 6
Get the sum of 2^(-i) with i from 0 to infinity:
In [3]: summation(2**-i, (i, 0, oo))
Out[3]: 2
Quickly compute the derivative of some hairy function:
In [5]: (x ** 2 * exp(x) * (sin(omega*x + phi))).diff(x)
Out[5]:
2 x 2 x x
omega*x *e *cos(omega*x + phi) + x *e *sin(omega*x + phi) + 2*x*e *sin(omega*x
+ phi)
If you have matplotlib installed then it can quickly plot things for
you. Otherwise it falls back to a great ascii poltting library:
In [11]: plot(x**2)
100 |
| . .
| . .
| . .
| . .
| . .
| . .
| .
50.0165 | --------.--------------------------------------..------
| . .
| . .
| . .
| .. ..
| . .
| .. ..
| .. ..
| .. ..
0.03305 | ..............
-10 0 10
Out[11]:
Plot object containing:
[0]: cartesian line: x**2 for x over (-10.0, 10.0)
Sympy is good for simple interactive work. It is also suitable for
complex problems but some of the interfaces are a little clunky when
you start needing to manipulate things programmatically. Here's a
non-trivial script that I used to compute the Butcher table for the
4th-order Gauss-Legendre method:
#!/usr/bin/env python
#
# Compute the Butcher table for the Gauss-Legendre 4th order scheme
from sympy import *
c1, c2, a1, a2, b1, b2, k1, k2, xn, h = symbols('c1 c2 a1 a2 b1 b2 k1 k2 xn h')
# Interpolation points
# Roots of the Legendre polynomial (transformed)
p1 = Rational(1, 2) - Rational(1, 6)*sqrt(3)
p2 = Rational(1, 2) + Rational(1, 6)*sqrt(3)
# Imagine that we are integrating the differential equation dxdt = x
# In that case the solution to the implicit problem for the stages is given by
# eq1 and eq2 below.
eq1 = h * (xn + a1*k1 +a2*k2) - k1
eq2 = h * (xn + b1*k1 +b2*k2) - k2
# Lets solve the implicit problem to get k1 and k2 in terms of the table
# coefficients.
soln = solve([eq1, eq2], [k1, k2])
# True solution of the differential equation
xnp1 = xn * exp(h)
# Estimate that we get from our integration step
xnp1bar = (xn + c1*soln[k1] + c2*soln[k2])
# The relative error between the true solution and estimate
E = ((xnp1bar - xnp1) / xnp1)
# Utility function to extract the terms from the Maclaurin series below
# Must be an easier way to do this.
def get_terms(expr, var):
def prep(expr):
expr = expr.subs(a2, p1 - a1).subs(b2, p2 - b1)
expr = collect(cancel(expand(expr)), h)
return expr
results = {}
order = -1
while order < 5:
coeff, order = expr.leadterm(h)
results[var**order] = prep(coeff)
expr = cancel(expr - coeff * var ** order)
return results
# We'll expand the error as a Maclaurin series and try to cancel the first
# three terms.
Eh = get_terms(series(E, h, n=6), h)
# Solve for the values of c, a1 and b1 that cancel the first three terms of
# the series expansion.
c_sol = solve([Eh[h**2], 1-c1-c2], [c1, c2])
for v in h**3, h**4:
Eh[v] = cancel(expand(Eh[v].subs(c1, c_sol[c1]).subs(c2, c_sol[c2])))
(a_sol, b_sol), = solve([Eh[h**3], Eh[h**4]], [a1, b1])
# Construct as a Matrix for pretty-printing
table = Matrix([
[p1 , a_sol, p1 - a_sol],
[p2 , b_sol, p2 - b_sol],
[nan, c_sol[c1], c_sol[c2]],
])
print('Butcher Table for Gauss-Legendre 4th order')
pprint(table)
The output (after ~5 seconds) is:
$ python gl4.py
Butcher Table for Gauss-Legendre 4th order
[ ___ ___ ]
[ \/ 3 1 \/ 3 1]
[- ----- + - 1/4 - ----- + -]
[ 6 2 6 4]
[ ]
[ ___ ___ ]
[ \/ 3 1 1 \/ 3 ]
[ ----- + - - + ----- 1/4 ]
[ 6 2 4 6 ]
[ ]
[ nan 1/2 1/2 ]
Compare this with the second table here:
http://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_method
Oscar
More information about the Tutor
mailing list