Discussion :

[prog_]

Bonjour,

je tente de résoudre un système d'équation différentielle de premier ordre sous Matlab et afficher l'évolution des variables en fonction du temps, mais j'obtient une matrice remplit de 'NaN' et le warning "Matrix is singular, close to singular or badly scaled.Results may be inaccurate. RCOND = NaN.'' Comment faire pour remedier à la situation ?
Voici mon code:

PFE.m:

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
 
global t;
global y;
% Conditions initiales
y0=zeros(47,1);
y0(1) = 0;  %EGFR_EGR
y0(2) = 0; %2EGFR_EGF
y0(3) = 0; %EGFR_P
y0(4) = 0; %Gr_G
y0(5) = 0; %PI3K_Gr_G
y0(6) = 0; %EGFR_PI3K; 
y0(7) = 0; %PIP3
y0(8) = 0; %PIP3_Akt_PDK
y0(9) = 0; %Akt_P
y0(10)= 0; %PIP2
y0(11)= 0; %EGFR_PL
y0(12)= 0; %EGFR_PLP
y0(13)= 0; %PLCr_P
y0(14)= 1.05*10^-7; %PLCr
y0(15)= 0; %PKC
y0(16)= 0; %Ca2+
y0(17)= 1.5*10^-7; %Shc
y0(18)= 0; %ShcP
y0(19)= 8.5*10^-7; %Grb2
y0(20)= 0; %PI3K*
y0(21)= 3.4*10^-8; %SOS
y0(22)= 0; %EGFR_G
y0(23)= 0; %EGFR_G_S
y0(24)= 0; %G_S
y0(25)= 0; %EGFR_Sh
y0(26)= 0; %EGFR_Sh_P
y0(27)= 0; %EGFR_Sh_G
y0(28)= 0; %EGFR_Sh_G_S
y0(29)= 0; %Sh_G
y0(30)= 0;%Sh_G_S
y0(31)= 1*10^-8; %PI3K
y0(32)= 1.107*10^-8; %Ras-GDP
y0(33)= 0 ;%Ras-GTP
y0(34)= 5.531*10^-9; %Raf
y0(35)= 0; %RafP
y0(36)= 8*10^-7; %PI
y0(37)= 1.9926*10^-7; %MEK
y0(38)= 0; %MEKP
y0(39)= 0; %MEKPP
y0(40)= 4.1514*10^-7; %ERK
y0(41)= 0; %ERKP
y0(42)= 0; %ERKPP
y0(43)= 3.33*10^-8; %EGF
y0(44)= 8.307*10^-9; %EGFR
y0(45)= 1*10^-8; %Akt
y0(46)= 4.3*10^-8; %Gab1
y0(47)= 1*10^-8; %Pkt
 
%Plage de temps de simulation 
tspan=[0, 100];
 
% Appelle au solveur ode
[t, y] = ode23s(@equation, tspan, y0);
 
plot(t,y(:,1));
plot(t,y(:,2));

et mon fichier equation.m contenant la fonction equation:

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
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Fonction equation.m
%Cette fonction contient l'ensemble des équations différentielles à
%résoudre, formé des réactions v1 à v44, déclaré ci-dessous.
 
function dydt=equation(t,y)
 
global t;
global y;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Constantes cinétiques
 
%
k1=3*10^6;k_1=0.06;
k2=1*10^7;k_2=0.1;
k3=1;k_3=0.01;
Vmax4=4.5*10^-7;Km4=5*10^-8;
k5=6*10^7;k_5=0.2;
k6=1;k_6=0.05;
k7=0.3;k_7=6*10^6;
Vmax8=4.5*10^-7;Km8=5*10^-8;
k9=3*10^6;k_9=0.05;
k10=1*10^7;k_10=0.06;
k11=0.03;k_11=6*10^7;
k12=0.015;k_12=1*10^5;
k13=9*10^7; k_13=0.6;
k14=6; k_14=0.06;
kcat15=0.2;Km15=8.3333*10^-7;
Vmax16=1.7*10^-9;Km16=3.4*10^-7;
k17=3*10^6; k_17=0.1;
k18=0.3;k_18=9*10^5;
k19=1*10^7; k_19=0.0214;
k20=0.12;k_20=2.4^5;
k21=3*10^6; k_21=0.1;
k22=3*10^6; k_22=0.064;
k23=01; k_23=2.1*10^7;
k24=9*10^6; k_24=0.0429;
kcat26=0.222 ;Km26=0.181;
Vmax28=0.289; Km28=0.0571;
k29=9.85 ;k_29=0.0985;
k30=3.6*10^7; k_30=0.05;
kcat31=1.53;Km31=11.7;
Vmax32=5*10^-6;Km32=3.3*10^-9;
kcat33=9*10^-9;Km33=4*10^-9;
Vmax35=5*10^-6;Km35=3.3*10^-7;
kcat37=0.138; Km37=6*10^-8;
Vmax39=6*10^-6;Km39=1*10^-7;
kcat41=16.67; Km41=5.54*10^-10;
kcat42=0.1667; Km42=5.4*10^-10;
k44=0.27; k_44=2;
kcat46=16.9; Km46=39.1;
k47=509; k_47=234;
Vmax48=2*10^4; kcat48=2*10^4; Km48=8*10^5;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Équations de réaction
%
 
v1=k1*y(43)*y(44)-k_1*y(1);
v2=k2*y(1)*y(1)-k_2*y(2);
v3=k3*y(2)- k_3*y(3);
v4=Vmax4*((y(3)/(Km4 * y(3))));
v5=k5*y(3)*y(14)-k_5*y(11);
v6=k6*y(11)-k_6*y(12);
v7=k7*y(12)-k_7*y(3)*y(13);
v8=Vmax8 *((y(13)/(Km8+y(13))));
v9=k9*y(3)*y(19)-k_9*y(22);
v10=k10*y(22)*y(21)-k_10*y(23);
v11=k11*y(23)-k_11*y(3)*y(24);
v12=k12*y(24)-k_12*y(19)*y(21);
v13=k13*y(3)-k_13*y(25);
v14=k14*y(25)-k_14*y(26);
v15=(kcat15*y(3)*y(17))/(Km15+y(17));
v16=Vmax16*(y(18)/(Km16+y(18)));
v17=k17*y(26)-k_17*y(27);
v18=k18*y(27)-k_18*y(3)*y(29);
v19=k19*y(27)*y(21)-k_19*y(28);
v20=k20*y(28)-k_20*y(3)*y(30);
v21=k21*y(18)*y(19)-k_21*y(29);
v22=k22*y(29)*y(21)-k_22*y(30);
v23=k23*y(30)-k_23*y(18)*y(24);
v24=k24*y(26)-k_24*y(28);
v26=(kcat26*y(28)*y(32))/(Km26+y(32));
v27=0;
v28=Vmax28*(y(33)/(Km28+y(33)));
v29=k29*y(6)-k_29*y(3)*y(20);
v31=(kcat31*y(33)*y(34))/(Km31+y(34));
v32=Vmax32*(y(35)/(Km32+y(35)));
v33=(kcat33*y(35)*y(37))/(Km33+y(37));
v34=(kcat33*y(38)*y(35))/(Km33+y(38));
v35=Vmax35*(y(38)/(Km35+y(38)));
v36=Vmax35*(y(39)/(Km35+y(39)));
v37=(kcat37*y(39)*y(40))/(Km33+y(40));
v38=(kcat37*y(39)*y(41))/(Km33+y(41));
v39=Vmax39*(y(41)/(Km39+y(41)));
v40=Vmax39*(y(42)/(Km39+y(42)));
v41=(kcat41*y(23)*y(32))/(Km41+y(32));
v43=0;
v44=k44*y(4)*y(31)-k_44*y(5);
v46=kcat46*y(20)*y(31)/Km46+y(20);
v47=k47*y(7)*y(45)*y(47)-k_47*y(8);
v48=Vmax48*y(8)/(Km48*(1+(y(9)/kcat48)+y(8)));
v49=0;
v50=0;
v51=0;
v52=v36;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Équations différentielles 
 
 
dydt=zeros(47,1);
 
dydt(1)=v1-v2;
dydt(2)=v2+v3-v4;
dydt(3)=v4-v3-v5-v9-v13-v27+v7+v11+v15+v18+v20+v29;
dydt(4)=v43-v44;
dydt(5)=v44-v27;
dydt(6)=v27-v46;
dydt(7)=v46-v47;
dydt(8)=v47-v8;
dydt(9)=v8;
dydt(10)=-v46-v49;
dydt(11)=v5-v6;
dydt(12)=v6-v7-v49;
dydt(13)=v7-v8;
dydt(14)=v8-v5;
dydt(15)=v49;
dydt(16)=v49;
dydt(17)=v16-v13;
dydt(18)=v15+v33-v16;
dydt(19)=v12-v9-v17-v21;
dydt(20)=v29-v46;
dydt(21)=v12-v10-v22;
dydt(22)=v9-v10;
dydt(23)=v10-v11-v41;
dydt(24)=v11-v12+v23;
dydt(25)=v13-v14;
dydt(26)=v14-v17-v15-v24;
dydt(27)=v17-v18-v19;
dydt(28)=v19+v24-v26-v20-v26;
dydt(29)=v21-v18-v22;
dydt(30)=v22-v20-v23;
dydt(31)=-v44;
dydt(32)=v26+v41-v28+v31;
dydt(33)=v28-v26;
dydt(34)=v31-v32;
dydt(35)=v31-v32-v33-v34;
dydt(36)=-v46;
dydt(37)=v35-v33;
dydt(38)=v33-v34-v35;
dydt(39)=v34-v36-v37-v38-v50;
dydt(40)=v39-v37;
dydt(41)=v37-v38-v39-v40;
dydt(42)=v38-v40;
dydt(43)=-v1;
dydt(44)=-v1;
dydt(45)=-v47;
dydt(46)=-v43;
dydt(47)=-v47;

Tout compile bien orsque je lance PFE.m, mais j'obtiens des NaN comme réponse ?

[salseropom]

Salut, je n'ai pas de solution précise à te proposer mais j'ai une idée. Quand tu veux résoudre un système linéaire Ax = b, et que ta matrice A est mal conditionnée, tu résouds en faite CAx = Cb ou C est une matrice "préconditionneur", de telle manière que la matrice CA soit "bien" conditionnée

Pour un système d'EDO, on pourrait reprendre la même idée :

x'(t) = f(t,x)

revient à résoudre

AX'(t) = A f(t,x)

A toi de trouver une "bonne" matrice A