How can I input sum over in function depend on time?
    3 views (last 30 days)
  
       Show older comments
    
I get error and the error is 
Index exceeds the number of array elements (1). I want to write this equation (blue highlight)

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 Comments
Accepted Answer
  Torsten
      
      
 on 14 Jul 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 Comments
More Answers (0)
See Also
Categories
				Find more on Parallel Computing Fundamentals in Help Center and File Exchange
			
	Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!








