Zane (and everyone else),
You might be interested in trying VFGEN (www.warrenweckesser.net/vfgen). I just released version 2.2.0. VFGEN is a program that takes an XML description of a vector field (i.e. differential equations) and generates code for a wide variety of numerical software packages, including SciPy.
A simple example is the damped pendulum:
è' = v
v' = -b*v/(m*L^2) - (g/L)*sin(è)
An XML vector field file for this system is
----------------------------------------------------------------
<?xml version="1.0" ?>
<VectorField
Name="pendulum"
Description="Pendulum Vector Field">
<Parameter
Name="g"
Description="gravitational constant"
DefaultValue="9.81" />
<Parameter
Name="b"
Description="friction constant"
DefaultValue="0.0" />
<Parameter
Name="L"
Description="pendulum length"
DefaultValue="1.0" />
<Parameter
Name="m"
Description="mass"
DefaultValue="1.0" />
<StateVariable
Name="theta"
Description="Angle, measured from straight down"
Formula="v"
PeriodFrom="0"
PeriodTo="2*Pi"
DefaultInitialCondition="Pi-0.01" />
<StateVariable
Name="v"
Description="angular velocity"
Formula="-b*v/(m*L^2)-(g/L)*sin(theta)"
DefaultInitialCondition="0.0" />
<Function
Name="energy"
Description="total energy (kinetic plus potential)"
Formula="m*L^2*v^2/2 - m*g*L*cos(theta)" />
</VectorField>
----------------------------------------------------------------
The VFGEN scipy command reads this file and generates this Python code:
----------------------------------------------------------------
#
# pendulum.py
#
# Python file for the vector field named: pendulum
# This file implements the vector field as
# the function pendulum_vf. This function
# can be used with the SciPy ODEINT function.
#
# This file was generated by the program VFGEN (Version:2.2.0)
# Generated on 17-Mar-2008 at 00:48
#
from math import *
import numpy
#
# The vector field.
#
def pendulum_vf(y_,t_,p_):
"""
The vector field function for the vector field "pendulum"
"""
Pi = pi
theta = y_[0]
v = y_[1]
g = p_[0]
b = p_[1]
L = p_[2]
m = p_[3]
f_ = numpy.zeros((2,))
f_[0] = v
f_[1] = -sin(theta)/L*g-1.0/m/(L*L)*b*v
return f_
#
# The Jacobian.
#
def pendulum_jac(y_, t_, p_):
"""
The Jacobian of the vector field "pendulum"
"""
Pi = pi
theta = y_[0]
v = y_[1]
g = p_[0]
b = p_[1]
L = p_[2]
m = p_[3]
# Create the Jacobian matrix:
jac_ = numpy.zeros((2,2))
jac_[0,1] = 1.0
jac_[1,0] = -1.0/L*cos(theta)*g
jac_[1,1] = -1.0/m/(L*L)*b
return jac_
#
# User function: energy
#
def pendulum_energy(y_, t_, p_):
"""
The user-defined function "energy" for the vector field "pendulum"
"""
Pi = pi
theta = y_[0]
v = y_[1]
g = p_[0]
b = p_[1]
L = p_[2]
m = p_[3]
return m*(L*L)*(v*v)/2.0-m*L*cos(theta)*g
----------------------------------------------------------------
VFGEN can also generate a simple command-line solver written in Python that uses the above function to solve the initial value problem. Take a look here: http://www.warrenweckesser.net/vfgen/menu_scipy.html
If you decide you would like to use a different package (e.g. MATLAB, Scilab, C with CVODE or with the GNU Scientific Library, Fortran with LSODA or RADAU5, etc), another one-line VFGEN command will generate the code that you need.
Best regards,
Warren Weckesser