[SciPy] Résolution d'un système d'équas diff

Version imprimable

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:

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 :

Citation:

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 !

17/07/2014, 22h29
Invité

Bonjour,

Merci d'utiliser le bouton '#' en haut à droite de l'éditeur de messages pour insérer du code. :tagcode:

@+.
20/07/2014, 14h00
captainspaulding

Excusez moi pour cette erreur de style mon bon seigneur. :aie:

Mon cas est-il si désespéré ?
20/07/2014, 14h38
Invité

Citation:

Envoyé par captainspaulding

Excusez moi pour cette erreur de style mon bon seigneur. :aie:

Ce n'est pas une erreur de style : Python est sensible aux indentations alors pour relire il vaut mieux avoir le code correctement formaté. ;)

Citation:

Mon cas est-il si désespéré ?

Absolument. :mrgreen:

Appelez le SAMU ! :mouarf:

@+.

Bon, je n'ai pas le background scientifique pour savoir ce qui bloque dans vos calculs, mais voici votre code remis en conformité avec les règles d'usage de Python :

Code:

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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
#!/usr/bin/env python
# -*- coding: utf-8 -*-
 
from math import exp
import numpy as np
from scipy.integrate import odeint
 
# main documentation
# PEP 8: http://legacy.python.org/dev/peps/pep-0008/
 
# Constantes
# cf http://legacy.python.org/dev/peps/pep-0008/#constants
 
K_GEL1 = 5.7e31
K_GEL2 = 3.1e14
K_GLC = 0.023
K_GLC2 = 2.9e-8
K_ALPHAMAL = 0.389
K_BETAMAL = 0.137
K_ALPHAMAL2 = 1.2e-7
K_BETAMAL2 = 8.4e-8
K_MLT = 0.117
K_MLT2 = 1.5e-8
K_DEX = 0.317
K_DEGALPHA = 6.9e30
E_DEGALPHA = 224.2
K_DEGBETA = 7.6e60
E_DEGBETA = 410.7
 
# fonctions
# cf http://legacy.python.org/dev/peps/pep-0008/#function-names
 
def temperature_profile (timeandtemp, t):
    """
        décrivez dans ce docstring ce que fait cette fonction
        e.g.
        ajustement du profil température en fct d'un abaque
        @timeandtemp et d'un instant @t;
    """
 
    # inits
 
    t = abs(float(t))
 
    tt00 = float(timeandtemp[0][0])
    tt01 = float(timeandtemp[0][1])
    tt10 = float(timeandtemp[1][0])
    tt11 = float(timeandtemp[1][1])
 
    if not t:
        T = 55.0
    elif t < 9.0:
        T = 55.0 + t
    elif t <= 9.0 + tt00:
        T = tt01
    elif t < 17.0 + tt00:
        T = tt01 + t
    elif t <= 17.0 + tt00 + tt01:
        T = tt11
    elif t < 23.0 + tt00 + tt10:
        T = tt11 + t
    else:
        T = float(timeandtemp[2][1])
    # end if
 
    return float(273.15 + T)
 
# end def
 
 
def flux (ci, temperature):
    """
        décrivez dans ce docstring ce que fait votre fonction
        e.g.
        calcul de profil de flux bla bla bla
    """
 
    # inits
 
    starchUG, starchG, glc, mal, mlt, dex, enzA, enzB = ci
 
    th = temperature
 
    # Relative activities
 
    if th < 336.15:
        a_alpha = -0.0011*th**3 + 1.091*th**2 - 352.89*th + 38008.3
        a_beta = 0.049*th - 13.9
    else:
        a_alpha = 0.0055*th**3 - 5.663*th**2 + 1941.9*th - 221864
        a_beta = 0.374*th + 128.3
    # end if
 
    # Equations
 
    if th < 333.15:
        r_gel = K_GEL1 * starchUG
    else:
        r_gel = K_GEL2 * starchUG
    # end if
 
    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 * exp(-E_DEGALPHA/(8.3145*th)) * enzA
    r_degB = K_DEGBETA * exp(-E_DEGBETA/(8.3145*th)) * enzB
    r_acA = K_DEGALPHA * exp(-E_DEGALPHA/(8.3145*th)) * a_alpha * enzA
    r_acB = K_DEGBETA * exp(-E_DEGBETA/(8.3145*th)) * 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]
    )
 
# end def
 
 
def sec_membre (ci, t, tempProf):
    """
        décrivez dans ce docstring ce que fait votre fonction
        e.g.
        calcul de bla bla bla
    """
 
    # inits
 
    temperature = temperature_profile(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
 
# end def
 
 
# Initial conditions
 
initCond = np.array([113.5, 0, 4, 5, 0, 0, 80000, 80000])
 
# Temperature steps in minutes
 
ts = temperature_steps = ((30, 64), (30, 72), (5, 78))
 
t0 = 0
 
tf = ts[0][0] + ts[1][0] + ts[2][0] + 23
 
nbPoints = 100
 
timeProfile = np.linspace(t0, tf, nbPoints)
 
print "\ntime profile:\n", timeProfile
 
print "\nrunning odeint():"
 
result = odeint(sec_membre, initCond, timeProfile, args=(temperature_steps,))
 
print "\nresult:\n", result

Peut-être qu'une fois simplifié de la sorte, vous verrez apparaître des erreurs que vous n'aviez pas vues auparavant ?

PEP 8 - bon usage de Python : http://legacy.python.org/dev/peps/pep-0008/

@+.

En optimisant les calculs dans def flux(), je me suis aperçu que la variable r_gel était définie mais jamais utilisée nulle part.

Quid est?

Code optimisé :

Code:

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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
#!/usr/bin/env python
# -*- coding: utf-8 -*-
 
from math import exp
import numpy as np
from scipy.integrate import odeint
 
# main documentation
# PEP 8: http://legacy.python.org/dev/peps/pep-0008/
 
# Constantes
# cf http://legacy.python.org/dev/peps/pep-0008/#constants
 
K_GEL1 = 5.7e31
K_GEL2 = 3.1e14
K_GLC = 0.023
K_GLC2 = 2.9e-8
K_ALPHAMAL = 0.389
K_BETAMAL = 0.137
K_ALPHAMAL2 = 1.2e-7
K_BETAMAL2 = 8.4e-8
K_MLT = 0.117
K_MLT2 = 1.5e-8
K_DEX = 0.317
K_DEGALPHA = 6.9e30
E_DEGALPHA = 224.2
K_DEGBETA = 7.6e60
E_DEGBETA = 410.7
 
# fonctions
# cf http://legacy.python.org/dev/peps/pep-0008/#function-names
 
def temperature_profile (timeandtemp, t):
    """
        décrivez dans ce docstring ce que fait cette fonction
        e.g.
        ajustement du profil température en fct d'un abaque
        @timeandtemp et d'un instant @t;
    """
 
    # inits
 
    t = abs(float(t))
 
    tt00 = float(timeandtemp[0][0])
    tt01 = float(timeandtemp[0][1])
    tt10 = float(timeandtemp[1][0])
    tt11 = float(timeandtemp[1][1])
 
    if not t:
        T = 55.0
    elif t < 9.0:
        T = 55.0 + t
    elif t <= 9.0 + tt00:
        T = tt01
    elif t < 17.0 + tt00:
        T = tt01 + t
    elif t <= 17.0 + tt00 + tt01:
        T = tt11
    elif t < 23.0 + tt00 + tt10:
        T = tt11 + t
    else:
        T = float(timeandtemp[2][1])
    # end if
 
    return float(273.15 + T)
 
# end def
 
 
def flux (ci, temperature):
    """
        décrivez dans ce docstring ce que fait votre fonction
        e.g.
        calcul de profil de flux bla bla bla
    """
 
    # inits
 
    starchUG, starchG, glc, mal, mlt, dex, enzA, enzB = ci
 
    th = temperature
 
    # Relative activities
 
    if th < 336.15:
        a_alpha = -0.0011*th**3 + 1.091*th**2 - 352.89*th + 38008.3
        a_beta = 0.049*th - 13.9
    else:
        a_alpha = 0.0055*th**3 - 5.663*th**2 + 1941.9*th - 221864.0
        a_beta = 0.374*th + 128.3
    # end if
 
    # Equations
 
    if th < 333.15:
        r_gel = K_GEL1 * starchUG   # r_gel defined but unused
    else:
        r_gel = K_GEL2 * starchUG   # r_gel defined but unused
    # end if
 
    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 + K_BETAMAL2*a_beta) * dex
    r_mlt2 = K_MLT2 * a_alpha * dex
    r_degA = K_DEGALPHA * exp(-E_DEGALPHA/(8.3145*th)) * enzA
    r_degB = K_DEGBETA * exp(-E_DEGBETA/(8.3145*th)) * enzB
    r_acA = r_degA * a_alpha
    r_acB = r_degB * a_beta
 
    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
    ])
 
# end def
 
 
def sec_membre (ci, t, tempProf):
    """
        décrivez dans ce docstring ce que fait votre fonction
        e.g.
        calcul de bla bla bla
    """
 
    # inits
 
    temperature = temperature_profile(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
 
# end def
 
 
# Initial conditions
 
initCond = np.array([113.5, 0, 4, 5, 0, 0, 80000, 80000])
 
# Temperature steps in minutes
 
ts = temperature_steps = ((30, 64), (30, 72), (5, 78))
 
t0 = 0
 
tf = ts[0][0] + ts[1][0] + ts[2][0] + 23
 
nbPoints = 100
 
timeProfile = np.linspace(t0, tf, nbPoints)
 
print "\ntime profile:\n", timeProfile
 
print "\nrunning odeint():"
 
result = odeint(sec_membre, initCond, timeProfile, args=(temperature_steps,))
 
print "\nresult:\n", result

@+.