Ode solver, how to pass V=0:0.1:1 used in the differential equations

1 vue (au cours des 30 derniers jours)
El Houssain Sabour
El Houssain Sabour le 20 Déc 2022
Commenté : Sam Chak le 20 Déc 2022
The system of differential equations of the following format:
function dy=ff_corr(t,y)
dy=zeros(7,1);
E_B1 = 0.7; alpha_B1= 0.5; E_B2 = 0.8; alpha_B2 = 0.5; E_A1 = 0.2;alpha_A1 = 0.5; E_A2 = 0.6;
alpha_A2 = 0.5; E_A3 = 0.9; alpha_A3 = 0.5;E_C7 = 0.65; alpha_C7 = 0.5;k_rev = 0.3; T = 343.15;
F = 96487;R = 8.31434;r_ox = 0;theta = 1;b = (R.*T)./F;TauPt = 2.15; TauC = 4.6;
dy(1) =(F./TauPt).*(4.7e-5.*(y(3).*exp(alpha_B1.*(V-E_B1-r_ox.*theta)./b)-(k_rev.*y(1).*exp(-(1- alpha_B1).*(V-E_B1)./b)))-1.6e-1.*(y(1).*exp(alpha_B1.*(V-E_B2-r_ox.*theta)./b)-(k_rev.*y(2).*exp(-(1-alpha_B2).*(V-E_B2)./b)))-9.3e-3.*y(1).*y(4).*exp(alpha_C7.*(V-E_C7)./b));
dy(2) =(F./TauPt).*(1.5e-1.*(y(1).*exp(alpha_B1.*(V-E_B2-r_ox.*theta)./b)-(k_rev.*y(2).*exp(-(1-alpha_B2).*(V-E_B2)./b))));
dy(3) = (F./TauPt).*(9.3e-3.*y(1).*y(4).*exp(alpha_C7.*(V-E_C7)./b)-4.7e-5.*(y(3).*exp(alpha_B1.*(V-E_B1-r_ox.*theta)./b)-(k_rev.*y(1).*exp(-(1-alpha_B1).*(V-E_B1)./b))));
dy(4) =(F./TauC).*(1.62e-8.*((y(6).*exp(alpha_A1.*((V-E_A1)./b)))-(y(4).*exp(-(1-alpha_A1).*((V-E_A1)./b))))-6.7e-1.*((y(4).*exp(alpha_A2.*((V-E_A2)./b)))-(y(5).*exp(-(1-alpha_A2).*((V-E_A2)./b))))-1.17e-1.*y(4).*exp(alpha_A3.*((V-E_A3)./b))-9.3e-3.*y(1).*y(4).*exp(alpha_C7.*(V-E_C7)./b));
dy(5) =(F./TauC).*6.6e-1.*((y(4).*exp(alpha_A2.*((V-E_A2)./b)))-(y(5).*exp(-(1-alpha_A2).*((V-E_A2)./b))));
dy(6) =(F./TauC).*(9.3e-3.*y(1).*y(4).*exp(alpha_C7.*(V-E_C7)./b)-1.62e-8.*(y(6).*exp(alpha_A1.*(V-E_A1)./b))-(y(4).*exp(-(1-alpha_A1).*(V-E_A1)./b))+1.17e-1.*y(4).*exp(alpha_A3.*((V-E_A3)./b)));
dy(7) =-3.14e-14.*y(3).*(exp(0.5.*(V-1.188)./0.3)-(5e-5.*exp(-(1-0.5).*(V-1.188)./0.3)));
end
y0=[0;0;1;0;0;1;2e-9];
V=0:0.1:10;
[t, y] = ode45(@ff_corr,[0 10], y0);
  1 commentaire
Sam Chak
Sam Chak le 20 Déc 2022
Does the system have an equilibrium when ? Just curious!
You forgot to supply the time vector for V. Is supposed to look like this? Something that adds "energy" into the system over time?
V = 0:0.1:10;
t = V;
plot(t, V, '.'), grid on, xlabel('t'), ylabel('V(t)')

Connectez-vous pour commenter.

Réponses (2)

VBBV
VBBV le 20 Déc 2022
This is one way
clearvars
y0=[0;0;1;0;0;1;1];
V=0:0.1:1;
for k = 1:length(V)
[t, y] = ode45(@(t,y) ff_corr(t,y,V(k)),[0 10], y0);
plot(t,y(:,1))
hold on
end
function dy=ff_corr(t,y,V)
dy=zeros(7,1);
E_B1 = 0.7; alpha_B1= 0.5; E_B2 = 0.8; alpha_B2 = 0.5; E_A1 = 0.2;alpha_A1 = 0.5; E_A2 = 0.6;
alpha_A2 = 0.5; E_A3 = 0.9; alpha_A3 = 0.5;E_C7 = 0.65; alpha_C7 = 0.5;k_rev = 0.3; T = 343.15;
F = 96487;R = 8.31434;r_ox = 0;theta = 1;b = (R.*T)./F;TauPt = 2.15; TauC = 4.6;
dy(1) =(F./TauPt).*(4.7e-5.*(y(3).*exp(alpha_B1.*(V-E_B1-r_ox.*theta)./b)-(k_rev.*y(1).*exp(-(1-alpha_B1).*(V-E_B1)./b)))-1.6e-1.*(y(1).*exp(alpha_B1.*(V-E_B2-r_ox.*theta)./b)-(k_rev.*y(2).*exp(-(1-alpha_B2).*(V-E_B2)./b)))-9.3e-3.*y(1).*y(4).*exp(alpha_C7.*(V-E_C7)./b));
dy(2) =(F./TauPt).*(1.5e-1.*(y(1).*exp(alpha_B1.*(V-E_B2-r_ox.*theta)./b)-(k_rev.*y(2).*exp(-(1-alpha_B2).*(V-E_B2)./b))));
dy(3) = (F./TauPt).*(9.3e-3.*y(1).*y(4).*exp(alpha_C7.*(V-E_C7)./b)-4.7e-5.*(y(3).*exp(alpha_B1.*(V-E_B1-r_ox.*theta)./b)-(k_rev.*y(1).*exp(-(1-alpha_B1).*(V-E_B1)./b))));
dy(4) =(F./TauC).*(1.62e-8.*((y(6).*exp(alpha_A1.*((V-E_A1)./b)))-(y(4).*exp(-(1-alpha_A1).*((V-E_A1)./b))))-6.7e-1.*((y(4).*exp(alpha_A2.*((V-E_A2)./b)))-(y(5).*exp(-(1-alpha_A2).*((V-E_A2)./b))))-1.17e-1.*y(4).*exp(alpha_A3.*((V-E_A3)./b))-9.3e-3.*y(1).*y(4).*exp(alpha_C7.*(V-E_C7)./b));
dy(5) =(F./TauC).*6.6e-1.*((y(4).*exp(alpha_A2.*((V-E_A2)./b)))-(y(5).*exp(-(1-alpha_A2).*((V-E_A2)./b))));
dy(6) =(F./TauC).*(9.3e-3.*y(1).*y(4).*exp(alpha_C7.*(V-E_C7)./b)-1.62e-8.*(y(6).*exp(alpha_A1.*(V-E_A1)./b))-(y(4).*exp(-(1-alpha_A1).*(V-E_A1)./b))+1.17e-1.*y(4).*exp(alpha_A3.*((V-E_A3)./b)));
dy(7) =-3.14e-14.*y(3).*(exp(0.5.*(V-1.188)./0.3)-(5e-5.*exp(-(1-0.5).*(V-1.188)./0.3)));
end
  2 commentaires
El Houssain Sabour
El Houssain Sabour le 20 Déc 2022
first of all thank you, unfortunately I got this error message.
Error using horzcat
Requested 6x178957200 (8.0GB) array exceeds maximum array size preference (8.0GB). This might cause MATLAB to become unresponsive.
Error in ode45 (line 476)
yout = [yout, zeros(neq,chunk,dataType)]; %#ok<AGROW>
Error in untitled4 (line 5)
[t, y] = ode45(@(t,y) ff_corr(t,y,V(k)),[0 10], y0);
Related documentation

Connectez-vous pour commenter.


Torsten
Torsten le 20 Déc 2022
So V should grow linearily from 0 to 10 in the differential equation ?
Then define it in @ff_corr as V = 10*t.

Catégories

En savoir plus sur Programming dans Help Center et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by