[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Calling lsode directly from a c++ program via liboctave?
From: |
Douglas Eck |
Subject: |
Re: Calling lsode directly from a c++ program via liboctave? |
Date: |
Tue, 23 Oct 2001 12:31:51 +0200 |
User-agent: |
Mozilla/5.0 (X11; U; Linux i686; en-US; rv:0.9.4) Gecko/20010913 |
I figured it out. I think it's kind of cool. And fast.
Attached is an example that shows an example
of calling the ODE solver lsode directly in
a dynamically-linked c++ file as opposed to calling
it from octave.
The example used is a stiff ODE in the form of a
Fitzhugh-Nagumo relaxation oscillator. This oscillator
is a model of a neuron being driven with constant voltage
from a "space clamp" on it's axon. It is a two-varible
simplification of the famous Hodgkin Huxley equations.
If anyone can improve on it, please let me know.
Cheers,
Doug
fitz_driver.m
Description: type
/*
To make, type mkoctfile fitz.cc
This was tested using octave 2.1.43. I have
no idea if it works in 2.0
The example used is a stiff ODE in the form of a
Fitzhugh-Nagumo relaxation oscillator. This oscillator
is a model of a neuron being driven with constant voltage
from a "space clamp" on it's axon. It is a two-varible
simplification of the famous Hodgkin Huxley equations.
The Fitzhugh-Nagumo equations were derived simultaneoulsly by Fitzhugh and
Nagumo et.al.
R. Fitzhugh, Impulses and physiological states in theoretical
models of nerve membranes, 1182-Biophys. J., 1 (1961), pp. 445--466. and
J. Nagumo, S. Arimoto, and S. Yoshizawa, An active pulse
transmission line simulating 1214-nerve axons, Proc. IRL, 50 (1960), pp.
2061--2070.
Douglas Eck address@hidden
This code borrows heavily from John Eaton's lsode.cc
*/
#include <octave/config.h>
#include <octave/LSODE.h>
#include <octave/variables.h>
#include <octave/oct.h>
#include <octave/ov-struct.h>
#define FITZ_THOLD .2
#define FITZ_CVOLT .112
#define FITZ_SHUNT 1.2
static LSODE_options lsode_opts;
static ColumnVector fn_ode_epsilons;
ColumnVector fn_ode_helper (const ColumnVector& x, double t)
{
//This is the function called multiple times by lsode.
//Each oscillator has a voltage v and recovery w
//thus x is of length i=oscCount*2 with x(i)=v, x(i+1)=w
//The parameter epsilon controls the rate of oscillaton
//it is stored in a static column vector fn_ode_epsilons
int l=x.length();
ColumnVector xdot(l);
xdot.fill(0.0);
for (int i=0;i<l;i=i+2) {
xdot(i) = -x(i)*(x(i)-FITZ_THOLD) * (x(i)-1) - x(i+1) + FITZ_CVOLT;
xdot(i+1) = fn_ode_epsilons(i/2) * (x(i)-FITZ_SHUNT * x(i+1));
}
return xdot;
}
DEFUN_DLD (fitz, args,,
"Fitzhugh-Nagumo oscillator solved using lsode.
Calling syntax: fitz(x0,epsilons,t) where x0
is a vector defining starting conditions
of the oscillators. Each oscillator is represented
using a duple (v,w) v=starting voltage
w=starting recovery, with initial conditions in
x0=(v1,w1,v2,w2,...,vK,wK) where K=oscCount
Epsilon is the variable controlling the rate (period) of the
oscillator. It is stored in a static ColumnVector.
")
{
octave_value_list retval;
if (args.length()>2 | args.length()<2) {
print_usage("fitz");
}
ColumnVector state (args(0).vector_value ());
fn_ode_epsilons=ColumnVector(args(1).vector_value());
ColumnVector out_times (args(2).vector_value ());
double tzero = out_times (0);
int nsteps = out_times.capacity ();
ODEFunc func (fn_ode_helper);
LSODE ode (state, tzero, func);
ode.copy (lsode_opts);
int nstates = state.capacity ();
Matrix output (nsteps, nstates + 1);
cout << "Computing Fitzhugh-Nagumo for times " << out_times(0)
<< " to " << (int)out_times(out_times.length()-1) << endl;
output = ode.integrate (out_times);
retval.resize (1);
retval(0) = output;
return retval;
}