Bonjour à tous,

je suis en train d'essayer de modéliser avec Python la dégradation enzymatique de l'amidon et je bloque à la résolution du système. Avec ce code :

Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
 
#!/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,))
j'obtiens cette erreur :

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.
J'avais l'habitude de résoudre ce type de problème sous matlab avec ode23tb mais avec Python je ne sais pas si ma méthode est bonne. Si quelqu'un est familier avec ce type de résolution votre aide me serait précieuse.

Merci !