Recently, I needed to make some plots with the solution of a system of ordinary differential equations (ODE) for a textbook. I know a great tool to easily make beautiful scientific plots — matplotlib/pylab (www.matplotlib.sf.net). A very simplistic function for integrating ODE comes with it — rk4. I know, that there are better tools for python provided with scipy, but it is not installed on my workstation (which I don’t have root access to), I am too lazy to compile it locally.
OK, so we have this rk4 function. It takes 3 arguments: a function dx/dt = f(x,t), x0 (x at time zero) and time points, t. In the case of a system of ODEs, y x, x0 and t are vectors of floats. First argument, a function, can take only two arguments: x and current time. Not that much flexibility, eh?
Naturally, one would want to set some parameters for the ODE system. How to capture/pass them? First, I saw the solution to make a class, which would store parameters and provide a method for integration:
## -- `Dumb class style' --
## Call example:
##
##>> fhn_slow_instance = FH_N_slow()
##>> rk4(fhn_slow_instance.system, state_zero, time_vector)
##
## where rk4 is a simple ODE integrator from pylab/matplotlib
class FH_N_slow:
def __init__(self, I = 0.6, p=(1.0, 1.7, 0.5)):
self.I = I
self.p = p
def dv(v,w): return v - v**3/3. - w + self.I
def dw(v,w): return 0.08*(polyval(self.p,v) - w)
def system(self, state, time):
v,w = state
return self.dv(v,w), self.dw(v,w)
FH_N_slow — is a modification of the FitzHugh-Nagumo model of a firing neuron which allows for long interspike intervals thanks to parabolic nullcline for w and a `narrow’ channel between the two nullclines. The original idea of that was coined by Hindmarsh and Rose some 20+ years before, but they had a strange form of the FitzHugh equations (not like the canonical representation).
This example seems to be OK, but can we simplify it? It seems natural to define a function within a closure of parameters to capture them. So, we need closures. I ended up with something like this:
## -- `Functional style' --
## Doesn't it look like a closure in Lisps?
## Call example:
##
##>> rk4(fhn_slow1(), state_zero, time_vector)
def fhn_slow1(I=0.6, p = (1.0, 1.7, 0.5)):
def dv(v,w): return v - v**3/3. - w + I
def dw(v,w): return 0.08*(polyval(p,v) - w)
def system(state, time):
v,w = state
return dv(v,w), dw(v,w)
return wrapper
Our function fhn_slow1 returns another function, closed over the values for I and p.
This looks neat to me.
And here is the phase portrait and the time realisation of the v (membrane potential) variable:
