Series expansion with swiginac (linux)

Gribouillis 1 Tallied Votes 479 Views Share

This script computes the formal series expansion of a mathematical function on the command line using the python module swiginac, an interface to the CAS Ginac. Typical invocation in a linux terminal looks like

$ serexp.py -n 7 "log(1+x)*(1+x+x**2)**(-1)" 
1*x+(-3/2)*x**2+5/6*x**3+5/12*x**4+(-21/20)*x**5+7/15*x**6+Order(x**7)

As far as I know, swiginac does not work in windows, so use a virtual (k)ubuntu machine in this case.

#!/usr/bin/env python
# -*-coding: utf8-*-
# DATE: 2014-May-13
# AUTHOR: Gribouillis for the python forum at www.daniweb.com
# LICENSE: Public Domain
# PYTHON: 2.7

from __future__ import (absolute_import, division,
                        print_function, unicode_literals)
__version__ = "0.1.0"
__doc__ = '''Script to compute a formal series expansion on the command line.
Invoke with option -h for help.

Tested in Kubuntu Linux, after installing swiginac with

    $ sudo aptitude install python-swiginac
'''
import argparse
import ast
import swiginac

class V(ast.NodeVisitor):
    def __init__(self):
        self.names = set()
    def visit_Name(self, node):
        self.names.add(node.id)

def my_eval(expr, namespace):
    assert 'x' in namespace
    node = ast.parse(expr)
    visitor = V()
    visitor.visit(node)
    names = visitor.names
    for name in names:
        if name == 'x':
            continue
        try:
            namespace[name] = getattr(swiginac, name)
        except AttributeError:
            raise RuntimeError(('Invalid name', name))
    return(eval(expr, namespace, namespace))

def main(args):
    inp = "exp( sqrt(1+sin(x)**3) ) - exp(1)"
    inp = args.expression
    point = args.x
    n = args.n
    D = dict()
    x = swiginac.symbol(b'x', b'x')
    D['x'] = x
    xarg = my_eval("x==({})".format(point), D)
    expr = my_eval(inp, D)
    print(expr.series(xarg, n))
    
if __name__ == '__main__':
    parser = argparse.ArgumentParser(description=
        """Compute a formal series expansion with swiginac.
        
for example:
        
$ PROG -n 5 "sqrt(1 + 2*x)"
1+1*x+(-1/2)*x**2+1/2*x**3+(-5/8)*x**4+Order(x**5)

Some expressions won't work, for example

"log(1+x)/(1+x+x**2)" raises an exception, while
"log(1+x)*(1+x+x**2)**(-1)" works as expected
""", formatter_class =argparse.RawTextHelpFormatter )
    parser.add_argument('expression',
        action = 'store',
        help = 'function to expand such as "sqrt(1+sin(Pi*x)**3)"',
        metavar = 'EXPRESSION',
    )
    parser.add_argument('-x',
        action = 'store',
        dest = 'x',
        help = 'point (default to 0). Expressions are allowed, such as "acos(0)"',
        metavar = 'POINT',
        default = '0',
    )
    parser.add_argument('-n',
        action = 'store',
        dest = 'n',
        help = 'order (defaults to 5). Integer.',
        default = 5,
        type = int,
        metavar = 'ORDER',
    )
    args = parser.parse_args()
    main(args)