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

1 回表示 (過去 30 日間)
El Houssain Sabour
El Houssain Sabour 2022 年 12 月 20 日
コメント済み: Sam Chak 2022 年 12 月 20 日
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 件のコメント
Sam Chak
Sam Chak 2022 年 12 月 20 日
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)')

サインインしてコメントする。

回答 (2 件)

VBBV
VBBV 2022 年 12 月 20 日
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 件のコメント
El Houssain Sabour
El Houssain Sabour 2022 年 12 月 20 日
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

サインインしてコメントする。


Torsten
Torsten 2022 年 12 月 20 日
So V should grow linearily from 0 to 10 in the differential equation ?
Then define it in @ff_corr as V = 10*t.

カテゴリ

Help Center および File ExchangeOrdinary Differential Equations についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by