用SciPy求解一组常微分方程

Question

我正在尝试解决一组常微分方程，以模拟淀粉在淀粉酶（酶的一种）作用下的水解过程。当我尝试解这组方程时，我遇到了一个

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可能不是合适的解法，但我不知道该用哪个。也许有人对这种类型的方程比较熟悉？

scipy 数值求解常微分方程 odeint 水解过程

用SciPy求解一组常微分方程

1 个回答

撰写回答