How to use 5 function coupled each other using ODE45? it is possible?
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
cindyawati cindyawati
el 20 de Jun. de 2023
Comentada: cindyawati cindyawati
el 20 de Jun. de 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 comentarios
Torsten
el 20 de Jun. de 2023
Editada: Torsten
el 20 de Jun. de 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
Respuesta aceptada
Alan Stevens
el 20 de Jun. de 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 comentarios
Alan Stevens
el 20 de Jun. de 2023
Yes, you can also use ode45 - it just uses more points over the time span.
Más respuestas (0)
Ver también
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!