【发布时间】：2014-09-07 15:35:08
【问题描述】：
我正在尝试求解一组 ODE，以模拟淀粉酶（酶）对淀粉的水解。当我尝试求解方程组时，我得到一个
lsoda-- at current t (=r1), mxstep (=i1) steps
taken on this call before reaching tout
In above message, I1 = 500
In above message, R1 = 0.6333483931400E+00
Excess work done on this call (perhaps wrong Dfun type).
Run with full_output = 1 to get quantitative information.
错误。代码：
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import numpy as np
import math
from scipy.integrate import odeint
def TemperatureProfile (timeandtemp,t) :
if t==0.0 :
T = 273.15+55
elif t>0.0 and t<9.0 :
T = 273.15+55+t
elif t>= 9.0 and t<=9.0+timeandtemp[0][0] :
T = 273.15+ timeandtemp[0][1]
elif t>9.0+timeandtemp[0][0] and t<9.0+timeandtemp[0][0]+8.0 :
T = 273.15+ timeandtemp[0][1]+t
elif t>=9.0+timeandtemp[0][0]+8.0 and t<=9.0+timeandtemp[0][0]+8.0+timeandtemp[1][0] :
T = 273.15+timeandtemp[1][1]
elif t>9.0+timeandtemp[0][0]+8.0+timeandtemp[1][0] and t<9.0+timeandtemp[0][0]+8.0+timeandtemp[1][0]+6:
T = 273.15+ timeandtemp[1][1]+t
elif t>=9.0+timeandtemp[0][0]+8.0+timeandtemp[1][0]+6 and t<9.0+timeandtemp[0][0]+8.0+timeandtemp[1][0]+6 :
T = 273.15+timeandtemp[2][1]
return T
def Flux (ci,temperature) :
'''
Constants
'''
K_gel1 = 5.7*10**31
K_gel2 = 3.1*10**14
k_glc = 0.023
k_glc2 = 2.9*10**-8
k_alphamal = 0.389
k_betamal = 0.137
k_alphamal2 = 1.2*10**-7
k_betamal2 = 8.4*10**-8
k_mlt = 0.117
k_mlt2 = 1.5*10**-8
k_dex = 0.317
k_degalpha = 6.9*10**30
E_degalpha = 224.2
k_degbeta = 7.6*10**60
E_degbeta = 410.7
starchUG = ci[0]
starchG = ci[1]
glc = ci[2]
mal = ci[3]
mlt = ci[4]
dex = ci[5]
enzA = ci[6]
enzB = ci[7]
'''
Relative activities
'''
if temperature < 336.15 :
a_alpha = -0.0011*temperature**3 + 1.091*temperature**2 - 352.89*temperature + 38008.3
a_beta = 0.049*temperature - 13.9
else :
a_alpha = 0.0055*temperature**3 - 5.663*temperature**2+ 1941.9*temperature- 221864
a_beta = 0.374*temperature + 128.3
'''
Equations
'''
if temperature < 333.15 :
r_gel = K_gel1*starchUG
else :
r_gel = K_gel2*starchUG
r_glc = k_glc*a_alpha*glc
r_mal = (k_alphamal*a_alpha+k_betamal*a_beta)*starchG
r_mlt = k_mlt*a_alpha*starchG
r_dex = k_dex*a_alpha*starchG
r_glc2 = k_glc2*a_alpha*dex
r_mal2 = k_alphamal2*a_alpha*dex+k_betamal2*a_beta*dex
r_mlt2 = k_mlt2*a_alpha*dex
r_degA = k_degalpha*math.exp(-E_degalpha/(8.3145*temperature))*enzA
r_degB = k_degbeta*math.exp(-E_degbeta/(8.3145*temperature))*enzB
r_acA = k_degalpha*math.exp(-E_degalpha/(8.3145*temperature))*a_alpha*enzA
r_acB = k_degbeta*math.exp(-E_degbeta/(8.3145*temperature))*a_beta*enzB
return np.array([r_glc,r_mal,r_mlt, r_dex, r_glc2, r_mal2, r_mlt2, r_degA, r_degB, r_acA, r_acB])
def Secmembre (ci,t,tempProf) :
temperature = TemperatureProfile(tempProf,t)
mS = np.mat([
[-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0],
[0, 0, 0, 0, 1, -1, -1, -1, 0, 0, 0]])
mS = np.transpose(mS)
mS = np.array(mS)
metabos = Flux(ci,temperature)
varMetabos = np.dot(metabos,mS)
return varMetabos
# Initial conditions
initCond = []
temperatureProfile = ((30,64),(30,72),(5,78)) # Temperature steps in minutes
t0 = 0
tf = temperatureProfile[0][0]+temperatureProfile[1][0]+temperatureProfile[2][0]+9+8+6
nbPoints = 100
timeProfile = np.linspace(t0,tf,nbPoints)
print timeProfile
initCond.append(113.5)
initCond.append(0)
initCond.append(4)
initCond.append(5)
initCond.append(0)
initCond.append(0)
initCond.append(80000)
initCond.append(80000)
initCond=np.array(initCond)
result = odeint(Secmembre,initCond,timeProfile,args=(temperatureProfile,))
在使用 Matlab 时，我曾经使用 ode23tb 求解。我认为 odeint 不是正确的求解器，但我不知道该使用哪一个。也许有人熟悉这种类型的方程？
【问题讨论】：
-
再次检查你的方程：如果你只解到
tf = 0.1，就没有错误。但是，解决方案呈指数增长，这不是物理的。您可能会收到错误，因为求解器发现解决方案最终会变得大于最大可能的浮点值（~ 1e308）。 -
特别是，
Secmembre返回 6 个变量的导数，但初始条件有 8 个变量。 （odeint 假设未指定的 2 个变量的导数为零）