How to use 5 function coupled each other using ODE45? it is possible?
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
cindyawati cindyawati
il 20 Giu 2023
Commentato: cindyawati cindyawati
il 20 Giu 2023
I have a coupled differential equation. I'm confused why my M3, O, and P values are 0. Even though the initial conditions M2, M3, O and P are 0 but only M2 has a value. Is there something wrong with the function?
function CM1 = mymode (t,M1,M2,M3,O,P)
M1= 10;
M2 = 0;
M3 = 0;
O = 0;
P=0;
delta=50;
gamma=75;
K1= 10^-4;
K2=5*10^-4;
K3=10^-3;
Ko=0.1;
n=3;
Oa=10;
Pa=100;
mu_1=10^-3;
mu_2=10^-3;
mu_3=10^-3;
mu_o=10^-4;
mu_p= 10^-5;
CM1= zeros(5,1);
CM1(1) = (delta*M1*(1-(M1/gamma))-2*K1*M1*M1-M1*(K2*M2)-((Oa-n)*K3*M1*M3)-((Pa-Oa)*Ko*M1*O)-(mu_1*M1));
CM1(2) = (K1*M1*M1)-(K2*M1*M2)-(mu_2*M2);
CM1(3) = (K2*M1*M2)-(K3*M1*M3)-(mu_3*M3);
CM1(4) = (K3*M1*M3)-(Ko*M1*O)-(mu_o*O);
CM1(5) = (Ko*M1*O)-(mu_p*P);
end
[t,M1,M2,M3,O,P] = ode45(@mymode, [0,100],[0,0.01])
plot (t,M1,M2,M3,O,P)
4 Commenti
Torsten
il 20 Giu 2023
Modificato: Torsten
il 20 Giu 2023
This is a Runge-Kutta-4 code for your problem. Try to understand how "runge_kutta_RK4" works on your system of equations to do better next time.
tstart = 0.0;
tend = 100.0;
dt = 0.01;
T = (tstart:dt:tend).';
Y0 = [10 0 0 0 0];
f = @myode;
Y = runge_kutta_RK4(f,T,Y0);
M1 = Y(:,1);
M2 = Y(:,2);
M3 = Y(:,3);
O = Y(:,4);
P = Y(:,5);
figure
subplot(3,1,1)
plot(T,M1),grid, title('M1')
subplot(3,1,2)
plot(T,M2),grid, title('M2')
subplot(3,1,3)
plot(T,M3),grid, title('M3')
figure
subplot(2,1,1)
plot(T,O),grid, title('O')
subplot(2,1,2)
plot(T,P),grid, title('P')
function Y = runge_kutta_RK4(f,T,Y0)
N = numel(T);
n = numel(Y0);
Y = zeros(N,n);
Y(1,:) = Y0;
for i = 2:N
t = T(i-1);
y = Y(i-1,:);
h = T(i) - T(i-1);
k0 = f(t,y);
k1 = f(t+0.5*h,y+k0*0.5*h);
k2 = f(t+0.5*h,y+k1*0.5*h);
k3 = f(t+h,y+k2*h);
Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
end
end
function CM1 = myode (~,MM)
M1 = MM(1);
M2 = MM(2);
M3 = MM(3);
O = MM(4);
P = MM(5);
delta=50;
gamma=75;
K1= 10^-4;
K2=5*10^-4;
K3=10^-3;
Ko=0.1;
n=3;
Oa=10;
Pa=100;
mu_1=10^-3;
mu_2=10^-3;
mu_3=10^-3;
mu_o=10^-4;
mu_p= 10^-5;
CM1= zeros(1,5);
CM1(1) = (delta*M1*(1-(M1/gamma))-2*K1*M1*M1-M1*(K2*M2)-((Oa-n)*K3*M1*M3)-((Pa-Oa)*Ko*M1*O)-(mu_1*M1));
CM1(2) = (K1*M1*M1)-(K2*M1*M2)-(mu_2*M2);
CM1(3) = (K2*M1*M2)-(K3*M1*M3)-(mu_3*M3);
CM1(4) = (K3*M1*M3)-(Ko*M1*O)-(mu_o*O);
CM1(5) = (Ko*M1*O)-(mu_p*P);
end
Risposta accettata
Alan Stevens
il 20 Giu 2023
Better like this:
MM0 = [10, 0, 0, 0, 0];
tspan = [0 100];
[t, MM] = ode15s(@mymode, tspan,MM0);
M1 = MM(:,1);
M2 = MM(:,2);
M3 = MM(:,3);
O = MM(:,4);
P = MM(:,5);
figure
subplot(3,1,1)
plot(t,M1),grid, title('M1')
subplot(3,1,2)
plot(t,M2),grid, title('M2')
subplot(3,1,3)
plot(t,M3),grid, title('M3')
figure
subplot(2,1,1)
plot(t,O),grid, title('O')
subplot(2,1,2)
plot(t,P),grid, title('P')
function CM1 = mymode (~,MM)
M1 = MM(1);
M2 = MM(2);
M3 = MM(3);
O = MM(4);
P = MM(5);
delta=50;
gamma=75;
K1= 10^-4;
K2=5*10^-4;
K3=10^-3;
Ko=0.1;
n=3;
Oa=10;
Pa=100;
mu_1=10^-3;
mu_2=10^-3;
mu_3=10^-3;
mu_o=10^-4;
mu_p= 10^-5;
CM1= zeros(5,1);
CM1(1) = (delta*M1*(1-(M1/gamma))-2*K1*M1*M1-M1*(K2*M2)-((Oa-n)*K3*M1*M3)-((Pa-Oa)*Ko*M1*O)-(mu_1*M1));
CM1(2) = (K1*M1*M1)-(K2*M1*M2)-(mu_2*M2);
CM1(3) = (K2*M1*M2)-(K3*M1*M3)-(mu_3*M3);
CM1(4) = (K3*M1*M3)-(Ko*M1*O)-(mu_o*O);
CM1(5) = (Ko*M1*O)-(mu_p*P);
end
3 Commenti
Alan Stevens
il 20 Giu 2023
Yes, you can also use ode45 - it just uses more points over the time span.
Più risposte (0)
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!