leakage in R-chain

From time to time I will post here some little code snippets from what I write.
All code is under GPL

This time it is a leakage resitance calculator:

(defun make-R-elem (R1 R2)
  (list R1 R2))

(defun axial-R (elem)
  (car elem))

(defun leak-R (elem)
  (cadr elem))

(defun make-R-chain (axial leak &key (number 10))
  (do ((chain nil)
       (n 0  (1+ n)))
      ((> n number) (nreverse chain))
    (push (make-R-elem axial leak) chain)))

(defun leak-resistance (chain)
 "Leakage resistanse of the repeated  circut:
       o---/\/\/\---+------o
         R1 (axial) |
                    /
                    \
                    / R2 (leakage)
                    \
                    /
                    |
                   --- Gnd
   elements are in a list: ((R1 R2) (R1 R2))
  "
  (let ((elem (car chain)))
   (if (null (cdr chain))
    (+ (axial-R elem) (leak-R elem))
    (+ (axial-R elem)
       (/ 1.0 (+ (/ 1.0 (leak-R elem))
                 (/ 1.0 (leak-resistance (cdr chain)))))))))

On how to define an ODE system in matplotlib

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: