matlab: comment résoudre les équations différentielles couplées dépendantes des données stockées dans les tableaux

Question

Je veux résoudre un système de deux équations différentielles ordinaires dans MATLAB. Les paramètres des ODE dépendent des données mesurées stockées dans deux tableaux F et T. Lorsque j'exécute le programme, j'obtiens toujours l'erreur indiquée ci-dessous. Je suis sûr que cela a quelque chose à voir avec les tableaux, parce que quand j'utilise des nombres uniques pour F et T (par exemple F = 60, T = 30;) le programme fonctionne bien.

Subscript indices must either be real positive integers or logicals. 

Error in dynamics (line 46) 
ddyn(1) = k1*F(t) + v_b(t) - k_1*dyn(1) - v_a(t); 

Error in ode23 (line 256) 
f(:,2) = feval(odeFcn,t+h*A(1),y+f*hB(:,1),odeArgs{:}); 

Error in main (line 33) 
[t,sol] = ode23(@dynamics , (1:1:3000),[0 0]); 

Voici le code que je l'utilise pour la fonction principale et le système ODE:

Fonction principale:

[t,sol] = ode45(@dynamics , (1:1:3000),[0 0]); 

système ODE:

function [ddyn] = dynamics(t,dyn) 

% constant numbers 
k1 = 10^-2; k_1 = 8* 10^-3; k2 = 10^-2; k_2 = 4*10^-3; 
V_max_a = 1.6; V_max_b = 3.5;  
K_M_a = 1.5*10^-3; K_M_b = 2*10^-3;  
K_a_F = 9.4*10^5; K_a_T = 3.9*10; K_b_F = 1.3*10^4; K_b_T = 1.2*10^-10; 
r_a_F = 4.3*10^7; r_a_T = 4.2*10^9; r_b_F = 6.9*10^-7; r_b_T = 6.2*10^-9;  

%arrays with data e.g. 
F = 1:3000; 
T = 1:3000; 

% programm works if I use numbers, e.g: 
%F = 60; T = 30; 


ddyn = zeros(2,1); 

R_a_F = (K_a_F + r_a_F* F)/(K_a_F + r_a_F); 
R_a_T = (K_a_T + r_a_T* T)/(K_a_T + r_a_T); 
R_b_F = (K_b_F + r_b_F* F)/(K_b_F + r_b_F); 
R_b_T = (K_b_T + r_b_T* T)/(K_b_T + r_b_T); 

v_a = (V_max_a*dyn(1))/(K_M_a + dyn(1))*R_a_F .*R_a_T; 
v_b = (V_max_b*dyn(2))/(K_M_b + dyn(2))*R_b_F .*R_b_T; 



ddyn(1) = k1*F(t) + v_b(t) - k_1*dyn(1) - v_a(t); 
ddyn(2) = k2*T(t) + v_a(t) - k_2*dyn(2) - v_b(t); 

Answer 1

Toutes les fonctions dans le Matlab ODE suite, incluant ode45, supposons têtre continu et utiliser un pas de temps dynamique pour atteindre un certain niveau de précision. En tant que tel, vous ne pouvez pas supposer t être un nombre entier et ne doit jamais être utilisé comme index comme vous le faites avec F(t). Pour citer the documentation:

If tspan contains more than two elements [t0,t1,t2,...,tf] , then the solver returns the solution evaluated at the given points. This does not affect the internal steps that the solver uses to traverse from tspan(1) to tspan(end) . Thus, the solver does not necessarily step precisely to each point specified in tspan .

Par conséquent, en supposant F et T sont des fonctions continues dans le temps, je vous recommande de faire une fonction qui effectue une interpolation dans le temps, plus que probable par interp1 et passer cette fonction à votre ODE fonctionne par parametrization. Par exemple:

tspan = 1:3000; 
Ffun = @(t) interp1(tspan,F,t); % Default is linear 
[t,sol] = ode45(@(t,dyn) dynamics(t,dyn,Ffun) , tspan , [0 0]); 

Ceci est juste un exemple mais devrait, espérons-le, être réparable.

---

En particulier, ode45 utilise la paire Dormand-Prince (4,5) Runkge-Kutta pour son intégration dans le temps; en bref, la fonction compare une solution de quatrième ordre et de cinquième ordre pour décider si le résultat de l'étape de temps actuelle est suffisant ou s'il devrait être réduit.