How can I input sum over in function depend on time?
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
cindyawati cindyawati
le 14 Juil 2023
Modifié(e) : Torsten
le 14 Juil 2023
I get error and the error is
Index exceeds the number of array elements (1). I want to write this equation (blue highlight)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1434418/image.png)
this is my code
tstart = 0;
tend = 100;
dt = 0.01;
T = (tstart:dt:tend).';
Y0 = [10 0 0 0 0 0 0 0];
f = @myode;
Y = fRK4(f,T,Y0);
M1 = Y(:,1);
M2 = Y(:,2);
M3 = Y(:,3);
M4 = Y(:,4);
M5 = Y(:,5);
M6 = Y(:,6);
O = Y(:,7);
P = Y(:,8);
figure
%subplot(3,1,1)
%hold on
plot(T,M1,'r','Linewidth',2)
xlabel('Times (days)')
ylabel('M1 (gr/ml)')
figure
%subplot(3,1,2)
%hold on
plot(T,M2,'b','Linewidth',2)
xlabel('Times (days)')
ylabel('M2 (gr/ml)')
figure
%subplot(3,1,3)
%hold on
plot(T,M3,'g','Linewidth',2)
xlabel('Times (days)')
ylabel('M3 (gr/ml)')
figure
%subplot(3,1,3)
%hold on
plot(T,M4,'b--','Linewidth',2)
xlabel('Times (days)')
ylabel('M4 (gr/ml)')
figure
%subplot(3,1,3)
%hold on
plot(T,M5,'r--','Linewidth',2)
xlabel('Times (days)')
ylabel('M5 (gr/ml)')
figure
%subplot(3,1,3)
%hold on
plot(T,M6,'g--','Linewidth',2)
xlabel('times (days)')
ylabel('M6 (gr/ml)')
figure
%subplot(2,1,1)
%hold on
plot(T,O,'k','Linewidth',2)
xlabel('Times (days)')
ylabel('O (gr/ml)')
figure
%subplot(2,1,2)
%hold on
plot(T,P,'m','Linewidth',2)
xlabel('Times (days)')
ylabel('P (gr/ml)')
function Y = fRK4(f,T,Y0)
N = numel(T);
n = numel(Y0);
Y = zeros(N,n);
Y(1,:) = Y0;
for j = 2:N
t = T(j-1);
y = Y(j-1,:);
h = T(j) - T(j-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(j,:) = y + h/6*(k0+2*k1+2*k2+k3);
end
K=2;
M=2;
sum= zeros(1,6);
for i = 2:6
summ=K(i).*M(i);
sum(i)=sum(i)+summ;
end
end
function CM1 = myode (~,MM)
M1 = MM(1);
M2 = MM(2);
M3 = MM(3);
M4 = MM(4);
M5 = MM(5);
M6 = MM(6);
O = MM(7);
P = MM(8);
delta=50;
gamma=75;
K1=10^-4;
K2=5*10^-4;
K3=10^-3;
K4=5*10^-3;
K5=10^-2;
K6=5*10^-2;
Ko=0.1;
n=6;
Oa=10;
Pa=100;
mu_1=10^-3;
mu_2=10^-3;
mu_3=10^-3;
mu_4=10^-3;
mu_5=10^-3;
mu_6=10^-3;
mu_o=10^-4;
mu_p= 10^-5;
CM1= zeros(1,8);
CM1(1) = (delta*M1*(1-(M1/gamma))-2*K1*M1*M1-M1*sum(i)-((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)-(K4*M1*M4)-(mu_4*M4);
CM1(5) = (K4*M1*M4)-(K5*M1*M5)-(mu_5*M5);
CM1(6) = (K5*M1*M5)-(K6*M1*M6)-(mu_6*M6);
CM1(7) = (K6*M1*M6)-(Ko*M1*O)-(mu_o*O);
CM1(8) = (Ko*M1*O)-(mu_p*P);
end
0 commentaires
Réponse acceptée
Torsten
le 14 Juil 2023
Never make changes in the RK4 function. Changes have to be made in the call to RK4 and in f.
tstart = 0;
tend = 100;
dt = 0.01;
T = (tstart:dt:tend).';
Y0 = [10 0 0 0 0 0 0 0];
f = @myode;
Y = fRK4(f,T,Y0);
M1 = Y(:,1);
M2 = Y(:,2);
M3 = Y(:,3);
M4 = Y(:,4);
M5 = Y(:,5);
M6 = Y(:,6);
O = Y(:,7);
P = Y(:,8);
figure
%subplot(3,1,1)
%hold on
plot(T,M1,'r','Linewidth',2)
xlabel('Times (days)')
ylabel('M1 (gr/ml)')
figure
%subplot(3,1,2)
%hold on
plot(T,M2,'b','Linewidth',2)
xlabel('Times (days)')
ylabel('M2 (gr/ml)')
figure
%subplot(3,1,3)
%hold on
plot(T,M3,'g','Linewidth',2)
xlabel('Times (days)')
ylabel('M3 (gr/ml)')
figure
%subplot(3,1,3)
%hold on
plot(T,M4,'b--','Linewidth',2)
xlabel('Times (days)')
ylabel('M4 (gr/ml)')
figure
%subplot(3,1,3)
%hold on
plot(T,M5,'r--','Linewidth',2)
xlabel('Times (days)')
ylabel('M5 (gr/ml)')
figure
%subplot(3,1,3)
%hold on
plot(T,M6,'g--','Linewidth',2)
xlabel('times (days)')
ylabel('M6 (gr/ml)')
figure
%subplot(2,1,1)
%hold on
plot(T,O,'k','Linewidth',2)
xlabel('Times (days)')
ylabel('O (gr/ml)')
figure
%subplot(2,1,2)
%hold on
plot(T,P,'m','Linewidth',2)
xlabel('Times (days)')
ylabel('P (gr/ml)')
function Y = fRK4(f,T,Y0)
N = numel(T);
n = numel(Y0);
Y = zeros(N,n);
Y(1,:) = Y0;
for j = 2:N
t = T(j-1);
y = Y(j-1,:);
h = T(j) - T(j-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(j,:) = y + h/6*(k0+2*k1+2*k2+k3);
end
end
function CM1 = myode (~,MM)
M1 = MM(1);
M2 = MM(2);
M3 = MM(3);
M4 = MM(4);
M5 = MM(5);
M6 = MM(6);
O = MM(7);
P = MM(8);
delta=50;
gamma=75;
K1=10^-4;
K2=5*10^-4;
K3=10^-3;
K4=5*10^-3;
K5=10^-2;
K6=5*10^-2;
Ko=0.1;
n=6;
Oa=10;
Pa=100;
mu_1=10^-3;
mu_2=10^-3;
mu_3=10^-3;
mu_4=10^-3;
mu_5=10^-3;
mu_6=10^-3;
mu_o=10^-4;
mu_p= 10^-5;
sumterm = K2*M2 + K3*M3 + K4*M4 + K5*M5;
CM1= zeros(1,8);
CM1(1) = (delta*M1*(1-(M1/gamma))-2*K1*M1*M1-M1*sumterm-((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)-(K4*M1*M4)-(mu_4*M4);
CM1(5) = (K4*M1*M4)-(K5*M1*M5)-(mu_5*M5);
CM1(6) = (K5*M1*M5)-(K6*M1*M6)-(mu_6*M6);
CM1(7) = (K6*M1*M6)-(Ko*M1*O)-(mu_o*O);
CM1(8) = (Ko*M1*O)-(mu_p*P);
end
2 commentaires
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Axis Labels dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!