[SciPy-User] Problem with combining Fsolve and Integrate.Odeint

[Johann Rohwer]

Sending again - my previous reply doesn't seem to have made it to the 
list...

On Tuesday 13 October 2009, Peter Cimermančič wrote:
> I'm trying to model system that is described with few ODEs.
> Function, where ODEs are in, is given as def function(y,t). It
> takes two arguments as you can see. y is an array of different
> species in the model, whereas t is an array of time steps. Then,
> I'd like to calculate steady state using fsolve, which takes
> function with one argument only. When trying to solve steady state,
> this error is raised: "TypeError: There is a mismatch between the
> input and output shape of diff_equations.". How could I solve my
> problem?

Here is a self-contained example of 2 simple ODEs with mass action 
kinetics, that calculates a time course with odeint and then uses the 
final point of the time course as an initial estimate for a steady-
state calculation with fsolve. As you can see, the arguments of the 
function de are handled OK between odeint and fsolve.

BTW if you are interested in this problem in a more general way, you 
might want to look at our PySCeS software (http://pysces.sf.net) which 
is for simulation of (bio)chemical kinetic networks. It has high-level 
functions to calculate time courses and steady states automatically 
(amongst others) and runs on top of scipy.

Regards
Johann

---------------------8-<--------------------------------

import scipy
scipy.pkgload('optimize', 'integrate')
import pylab

k1 = 10
k2 = 5
k3 = 8

def de(X,t):
    Xdot = scipy.zeros((2),'d')
    Xdot[0] = k1 - k2*X[0]
    Xdot[1] = k2*X[0] - k3*X[1]
    return Xdot

init_val = scipy.array([10.,0.1])
t_range = scipy.linspace(0,1,21)

t_course = scipy.integrate.odeint(de, init_val, t_range)
fin_t_course = scipy.copy(t_course[-1])

ss = scipy.optimize.fsolve(de, fin_t_course, args=None)
print ss

pylab.plot(t_range, t_course[:,0], t_range, t_course[:,1])
pylab.show()