News

Note: An error was discovered in the root handling. This is fixed in the latest version.

The wrapper is based on code posten on the cython-dev mailing list by Mr. Jon Olav Vik.

Highlights

• CVODE - Solver for stiff and nonstiff ordinary differential equation
• IDA - Solver for the solution of differential-algebraic equation (DAE) systems
• KINSOL - solver for nonlinear algebraic systems based on Newton-Krylov solver technology

The CVODE and IDA solvers support root finding and the solver throws an exception on finding a root. There is also an example of implementing the explicit equation in Cython for speed.

Bouncing ball example from OpenModelica manual

import array
import numpy as np
from sundials import *

class Ball:
    TINY = 1e-16
    def __init__(self, e = 0.7, g = 9.81):
        self.e = e
        self.g = g

    def f(self, t, y, sw):
        if sw[0]:
            return [0., 0.]
        else:
            return [-self.g, y[0]]
        
    def rootf(self, t, y, sw):
        return [y[1] + Ball.TINY]
    
    def solve(self, plot = True):
        t0 = 0
        y0 = array.array('d', [0.,1.])
        solver = CVodeSolver(RHS = self.f, ROOT = self.rootf, SW = [False],
                       abstol = 1.0e-6, reltol = 1.0e-6)

        solver.init(t0,y0)

        dt = .01
        iter = solver.iter(t0, dt)
        tres = array.array('d', [t0])
        hres = array.array('d', [y0[1]])
        while True:
            try:
                t, y = next(iter)
                tres.append(t)
                hres.append(y[1])
                    
            except CVodeRootException, info:
                if abs(info.y[0]) < 0.01:
                    solver.SW[0] = True
                    
                tres.append(info.t)
                hres.append(0.)
                
                solver.init(info.t, [-self.e*info.y[0], 0.])
            
            if t > 3.0:
                break
        
        if plot:
            import matplotlib
            import pylab

            pylab.plot(tres,hres, linewidth=1.0)
            pylab.xlabel('time (s)')
            pylab.ylabel('height (m)')
            pylab.title('Bouncing ball demo')
            pylab.grid(True)
            pylab.show()
        else:
            return tres, hres

ball = Ball()
ball.solve()