Bonjour, je cherche à faire coller un modèle biologiques à des valeurs expérimentales.
J'ai adapté un code existant à mon problème, mais visiblement ça ne marche pas.
J'obtiens le code erreur 21, index invalide.
Je ne vois pas quelle matrice je n'ai pas bien défini.
Code:
// Experimental data
clc
        // data from 25-02-13
y_exp(:,1) = [8000 23000 61000 101000 116000 192000 183000 196000]';
y_exp(:,2) = [13000 9000 53000 77000 116000 205000 232000 261000]';
y_exp(:,3) = [6000 26000 46000 67000 136000 192000 236000 303000]';
y_exp(:,4) = [8000 19000 8000 11000 8000 9000 17000 20000]';
y_exp(:,5) = [7000 16000 8000 20000 33000 40000 32000 85000]';
y_exp(:,6) = [9000 16000 15000 16000 36000 74000 57000 79000]';
        // data from 12-03-13
y_exp(:,7) = [11000 9000 51000 104000 128000 182000 192000 284000]';
y_exp(:,8) = [9000 7000 14000 26000 46000 100000 114000 162000]';
y_exp(:,9) = [11000 13000 74000 185000 226000 356000 228000 428000]';
y_exp(:,10) = [11000 24000 22000 30000 78000 112000 168000 156000]';
y_exp(:,11) = [13000 27000 111000 201000 222000 272000 316000 332000]';
y_exp(:,12) = [14000 25000 150000 260000 254000 364000 396000 452000]';
        // data from 27-03-13
y_exp(:,13) = [10000 10000 16000 33000 46000 84000 75000 104000]';
y_exp(:,14) = [7000 8000 40000 66000 135000 228000 180000 204000]';
y_exp(:,15) = [8000 12000 38000 78000 129000 180000 162000 248000]';
y_exp(:,16) = [6000 10000 11000 14000 11000 28000 23000 27000]';
y_exp(:,17) = [6000 7000 9000 25000 23000 31000 50000 49000]';
y_exp(:,18) = [5000 7000 22000 32000 42000 50000 55000 65000]';
        // data from 11-04-13
y_exp(:,19) = [7000 24000 116000 188000 224000 332000 328000 392857]';
y_exp(:,20) = [8000 15000 39000 118000 140000 196000 288000 427429]';
y_exp(:,21) = [8000 12000 43000 103000 144000 192000 300000 480857]';
y_exp(:,22) = [7000 9000 30000 39000 82000 128000 160000 232571]';
y_exp(:,23) = [5000 9000 25000 51000 75000 136000 196000 238857]';
y_exp(:,24) = [8000 11000 27000 45000 99000 148000 184000 257714]';

Fluxi=[148.27 148.71 140.59 51.65 54.36 52.76 148.27 148.71 140.59 233.60 238.9 244.0 148.27 148.71 140.59 51.65 47.12 43.29 241.30 241.80 242.50 ...
141.93 160.74 159.9]

DOi = [0.425 0.416 0.414 1.286 0.374 0.191 1.325 0.399 0.149 1.358 0.363 0.221 0.384 0.893 0.379 0.194 0.917 0.382 0.855 0.25 0.892 0.885 0.87 0.887]

// Model
function [dy] = mymodel(t,y,flag,u,Ks,DOi,Fluxi)
dy(1) = u*y(1)*y(2)/(y(2)+Ks)*t 
dy(2) = =Fluxi*(0.902623-0.02971*y(3)+0.0005421*y(3)*y(3)-0.00000557133*y(3)*y(3)*y(3)+0.0000000299771*y(3)*y(3)*y(3)*y(3) ...
-0.0000000000653467*y(3)*y(3)*y(3)*y(3)*y(3))
dy(3) = =12,789+54,7014*DOi+0,000196202*y(1)
dy = dy'; //transposition
endfunction
// starting guess
u = 0.25;
Ks = 50;
k = [u;Ks];
// --------------------------------------------------------------------------
function [val] = L_Squares(k)
val=[];
u = k(1);
Ks = k(2);
time = [0 2 4 6 8 10 12 14]';
y0 = y_exp(1,:);
t = time;
y_calc=ode(y0',0,time,mymodel);
y_calc = y_calc';
resid = (y_calc - y_exp).*(y_calc - y_exp);
val = sum(sum(resid));
endfunction
//---------------------------------------------------------------------------
tic
options = optimset("TolFun",1D-02,"TolX",1.D-02,"MaxFunEvals",1.0D+06,"MaxIter",1.0D+06);
[param,L_Squares,flag,details] = fminsearch(L_Squares,k,options);
elapsed_time = toc();
parametros = param'
u = parametros(1); Ks = parametros(2);
clc
if flag == 1 then
printf('Le système converge!!!\n\n');
else
printf('Le système ne converge pas!!!\n\n');
end;
printf('nombre d itération ---------------- %g\n',details.iterations);
printf('Somme des carrés ----------------- %g \n',L_Squares);
printf('Nombre d appels du modèle ------- %g\n',details.funcCount);
printf('Algorithme de minimisation------------ %s\n',details.algorithm);
printf('µ ');
write(%io(2),u)
printf('Ks :');
write(%io(2),Ks)
if elapsed_time<60 then
printf('Temps écoulé -------------------- %g s \n\n',elapsed_time);
else
printf('Temps écoulé -------------------- %g min \n\n',elapsed_time);
end;
time = [0 2 4 6 8 10 12 14]';
t = time;
y=ode([-1 1]',0,time,mymodel);
y = y';
clf
plot(time,y_exp(:,1),"ok")
set(gca(),"auto_clear","off")
plot(time,y_exp(:,2),"ob")
plot(time,y(:,1),"-k")
plot(time,y(:,2),"-Ks")
h1=legend('exp_data_1','exp_data_2','model','model')
Une des principales différences avec le code sur lequel je me suis appuyé est que j'ai 3 équations et 24 jeux de données à faire coller dessus contrairement à l'exemple avec 2 equations et 2 jeux de données.

Si vous pouviez éclairer ma lanterne, je vous en serais grandement reconnaissant!!
Merci