How to use 5 function coupled each other using ODE45? it is possible?
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
cindyawati cindyawati
am 20 Jun. 2023
Kommentiert: cindyawati cindyawati
am 20 Jun. 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 Kommentare
Torsten
am 20 Jun. 2023
Bearbeitet: Torsten
am 20 Jun. 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
Akzeptierte Antwort
Alan Stevens
am 20 Jun. 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 Kommentare
Alan Stevens
am 20 Jun. 2023
Yes, you can also use ode45 - it just uses more points over the time span.
Weitere Antworten (0)
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!