How to use 5 function coupled each other using ODE45? it is possible?
조회 수: 2 (최근 30일)
이전 댓글 표시
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
Torsten
2023년 6월 20일
편집: Torsten
2023년 6월 20일
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
채택된 답변
Alan Stevens
2023년 6월 20일
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
Alan Stevens
2023년 6월 20일
Yes, you can also use ode45 - it just uses more points over the time span.
추가 답변 (0개)
참고 항목
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!