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

[Warren Weckesser]

Hi Johann,

You can do what you want by using the `args` argument of both fsolve and 
odeint.

I have rewritten your example to do use `args`, and made several other 
changes that are more about my preferences than about getting it 
working.  In particular, if the differential equations depend on the 
parameters k1, k2, and k3, these should really be arguments to the 
function, not global variables.  Like fsolve, odeint also has an `args` 
argument that allows you to do this.

Heres the code:

----------

from numpy import empty, array, linspace
from scipy.optimize import fsolve
from scipy.integrate import odeint
import pylab


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


if __name__ == "__main__":
    k1 = 10
    k2 = 5
    k3 = 8

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

    t_course = odeint(de, init_val, t_range, args=(k1,k2,k3))
    fin_t_course = t_course[-1]

    ss = fsolve(de, fin_t_course, args=(0, k1, k2, k3))
    print ss

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

----------

Warren

Johann Rohwer wrote:
> 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()
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>   

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20091013/3b226d6b/attachment.html>