Swiginac - Extending Python with Symbolic Mathematics

Ola Skavhaug [1][2] Ondrej Certic [3][4]

Excerpt from the slides of a presentation of swiginac given at the Fenics'05 meeting in Chicago

GiNaC

GiNaC is not a CAS. GiNaC is a C++ library for applications in need of symbolic manipulation. Python is such an application.

Features:

Swiginac

Swiginac is a Python interface for GiNaC and can be imported as module:

>>> from swiginac import *

Symbols

Symbols are basic units in Swiginac:

>>> a = symbol('a', r'\alpha')
>>> b = symbol('beta')
>>> print a, b
a beta
>>> u = b + a
>>> print u
beta+a
>>> u.set_print_context('tex')
>>> print u
\beta+\alpha

All expressions in GiNaC are built with symbols. The drawback of this approach is that the level of abstraction is limited

Numbers

The numeric class can manipulate arbitrary precision integers in a very fast way. Rational numbers are automatically converted to fractions of coprime integers.

When specifying rational numbers, the expression that is the argument to numeric is evaluated using Pythons "standard" arithmetic, which is probabely not what we want:

>>> numeric(1/3)
0
>>> numeric(1./3)
0.33333333333333331483

To prevent evaluation by the "normal" Python operators, rational numbers can be input as string value or as (numerator, denominator):

>>> numeric("1/3")
1/3
>>> numeric(1, 3)
1/3
>>> numeric(3, 9)
1/3

Often, it is sufficient to specify one number in an expression as numeric:

>>> numeric(1)/3
1/3

Functions

Lots of functions are available

>>> sin(exp(b))
sin(exp(beta))
>>> v = tgamma(a+sqrt(b))
>>> print v.printlatex()
\Gamma(\alpha+\sqrt{\beta})

All trigonometric and hyperbolic functions are implemented in GiNaC, most of them interfaced in swiginac

Symbolic differentiation

Objects have the method diff for differentiation

>>> x = symbol('x'); y = symbol('y')
>>> P = x**5 + x**2 + y
>>> P.diff(x, 1)
2*x+5*x**4
>>> P.diff(x, 2)
2+20*x**3

diff exists as a function call too

>>> u = sin(exp(x))
>>> diff(u, x)
exp(x)*cos(exp(x))
>>> diff(u, x, 2)
exp(x)*cos(exp(x))-sin(exp(x))*exp(x)**2

Matrices

>>> mat1 = matrix(2,2)   # Two by two matrix
>>> mat1[0,0] = x
>>> mat1[1,1] = y
>>> mat1
[[x,0],[0,y]]
>>> mat1a = diag_matrix([x,y]) # Alternative definition
>>> print mat1.printlatex()
\left(\begin{array}{cc}x&0\\0&y\end{array}\right)
>>> mat2 = matrix([[sqrt(a),0],[1.0, cosh(b)]])
>>> print mat2.printc()
[[pow(a,(1.0/2.0)),0.0],[1.0000000000000000e+00,cosh(beta)]]

Simple integral support

We can construct integral objects and integrate either symbolically or numerically:

>>> integ = integral(x, 0, 1, x*x)
>>> integ
[integral object]
>>> print integ.printlatex()
\int_{0}^{1} dx\,x^{2}
>>> integ.eval_integ()
1/3
>>> integ.evalf()
0.33333333333333333332

Substitution

Algebraic objects in expressions can be substituted:

>>> u = sin(exp(b))
>>> v = u.subs(exp(b)==sqrt(a)) # v = sin(a**(1/2))
>>> w = v.subs(a==2).evalf()    # Convert sin(2**(1/2)) to `numeric`
>>> print v, w
sin(a**(1/2)) 0.98776594599273552706
>>> float(w)                    # Convert to Python double
0.98776594599273548

Sub-expressions do not match:

>>> x = symbol('x'); y = symbol('y'); z = symbol('z')
>>> u = sin(x+y+z)
>>> u.subs(x+y==4)
sin(z+y+x)
>>> u.subs([x==1, y==3])
sin(4+z)
>>> u.subs([x==1, y==2, z==3]) # Same as u.subs(x+y+z==6)
sin(6)

Solving linear systems

lsolve solves linear systems:

>>> x = symbol('x'); y = symbol('y')
>>> solution = lsolve([3*x + 5*y == 2, 5*x+y == -3], [x,y])
>>> solution
[x==-17/22, y==19/22]

Taylor series expansion

Expressions can expand themselves as a Taylor series:

>>> sin(x).series(x==0, 8)
1*x+(-1/6)*x**3+1/120*x**5+(-1/5040)*x**7+Order(x**8)
[1]Simula Research Laboratory
[2]Dept. of Informatics, University of Oslo
[3]Faculty of Mathematics and Physics
[4]Charles University in Prague