cython odeint example

In response to the request by Ryan I now post a full example of using cython and odeint together. This is a specially written example because the original work contains too much code from my work that I do not want to post on the web. The equations used in the example have been chosen rather randomly.

In this example I have created a sample Model class which is a callable class (you can use an instance of the class as the function parameter for odeint). The Model class is defined in model_cython.pyx and contains all the Cython code. In ode_script.py the model function is integrated with odeint and thus shows how to call the Cython class from a normal python script. Setup.py is needed to compile model_cython.pyx (see this post).

model_cython.pyx:

import sys

#define some constants a c floats but still callable from python scripts
#not really used here but left in here as an example
cpdef float R_GAS,K_B
R_GAS=8.31      #Gas constant J/molK
K_B=1.38066e-23 #Boltzmann constant J/K

#import some often used math functions from the c math library
cdef extern from "math.h":
    double sin(double)

cdef extern from "math.h":
    double fabs(double)

cdef extern from "math.h":
    double exp(double)

#define a function callable from c and python
cpdef float Velocity(float t,float amplitude=1.0,float period=4.):
    cdef f1
    f1=fabs(sin(t*period))*amplitude
    return(f1)

class Model:
    """
    This class can be used to integrate a differential equation
    """
    def __init__(self,float a,float p):
        """
        -a is the amplitude of the sine function that describes the 
            speed oscillation as a function time
        -p controls the period of the sine function
        """
        self.a=a
        self.p=p
        
    
    def __call__(self,float y,float t):
        """
        y is the current state of the initial value problem:
        y: position
        t: time
        """
        cdef float pos
        pos=y
        #Workhardening
        v=Velocity(t,amplitude=self.a,period=self.p)*exp(-pos)
        return(v)

ode_script.py:

#!/bin/env/python

from numpy import r_
import pylab,scipy.integrate,sys
sys.path.append(".")
from model_cython import *


#ampltiude and period
amp=1.0
period=4.0
#set initial position
pos0=0.0
#create instance of Model object with speed amplitude of 1
myModel=Model(amp,period)
#create time vector (4 seconds in 500 steps)
tvec=r_[0.0:4.0:500j]
#integrate ov
y=(scipy.integrate.odeint(myModel,pos0,tvec))
    

#plot results
pylab.figure(1)
pylab.clf()
pylab.subplot(211)
pylab.plot(tvec,y)
pylab.xlabel("Time (s)")
pylab.ylabel("Position (m)")
pylab.grid(b=True)
pylab.subplot(212)
v_list=[]
for t in tvec:
    v_list.append(Velocity(t,amplitude=amp,period=period))
pylab.plot(tvec,v_list)
pylab.xlabel("Time (s)")
pylab.ylabel("Undamped velocity (m/s)")

pylab.grid(b=True)
pylab.show()

setup.py:

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

ext_modules = [Extension("model_cython",["model_cython.pyx"],\
    libraries=["m"])]

setup(
    name= 'Generic model class',
    cmdclass = {'build_ext': build_ext},
    ext_modules = ext_modules
)

Hope this helps.