[SciPy-User] Problem with combining Fsolve and Integrate.Odeint
Warren Weckesser
warren.weckesser@enthought....
Tue Oct 13 10:44:42 CDT 2009
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@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
More information about the SciPy-User
mailing list