in this I apply a loop so that I will have different value of Y(:,11) for the different value of initial condition for Y(:,11), what is correct algorithm for this??

please check this code is it right or worng? and please lemme know he correct program for this.

ti = 0;
tf = 1E-3;
tspan=[ti tf];

k = 1:5000:2*pi;
for i = 1:length(k)
y0(i)=[1; 1; 1; 1; 1; 1; 1; 1; 1; 1; k(i); 0; 0; 0; 0];
N = 5;
tp = 5.4E-9;
[T,Y]= ode45(@(t,y) rate_eq(t,y,N),tspan,y0(i));
E(i) = mean(Y(:,11));
q(i) = (5.*(log(E(i))))./(2*pi);
end

Unable to perform assignment because the indices on the left side are not compatible with the size of the right side.

q

function dy = rate_eq(t,y,N)

dy = zeros(3*N,1);

P = 1.67;
a = 0.1;
tc = 230E-6;
tp = 5.4E-9;
o1 = 1E4;
o2 = 2E4;
o3 = 3E4;
o4 = 4E4;
o5 = 11E4;
k1 = 5.4E-3;
k2 = 10.8E-3;
k3 = 16.2E-3;
k4 = 22E-3;
k5 = 54E-3;

dy(1) = (P - y(1).*((abs(y(6)))^2 +1))./tc; 
dy(2) = (P - y(2).*((abs(y(7)))^2 +1))./tc; 
dy(3) = (P - y(3).*((abs(y(8)))^2 +1))./tc; 
dy(4) = (P - y(4).*((abs(y(9)))^2 +1))./tc; 
dy(5) = (P - y(5).*((abs(y(10)))^2 +1))./tc; 

dy(6)= (y(1)-a).*((y(6))./tp) + (k1./tp).*(y(7)).*cos(y(11)) + (k5./tp).*(y(10))*cos(y(15));
dy(7)= (y(2)-a).*((y(7))./tp) + (k2./tp).*(y(8)).*cos(y(12)) + (k1./tp).*(y(6))*cos(y(11));
dy(8)= (y(3)-a).*((y(8))./tp) + (k3./tp).*(y(9)).*cos(y(13)) + (k2./tp).*(y(7))*cos(y(12));
dy(9)= (y(4)-a).*((y(9))./tp) + (k4./tp).*(y(10)).*cos(y(14)) + (k3./tp).*(y(8))*cos(y(13));
dy(10)= (y(5)-a).*((y(10))./tp) + (k5./tp).*(y(9)).*cos(y(14)) +(k4./tp).* (y(6))*cos(y(15));

dy(11) = o1 - (k1./tp).*(y(6)/y(7)).*(sin(y(11))) - (k1./tp).*(y(10)/y(6)).*(sin(y(11))) + (k2./tp).*(y(8)./y(7))*sin(y(12)) + (k5./tp).*(y(10)/y(6)).*sin(y(15));
dy(12) = o2 - (k2./tp).*(y(7)/y(8)).*(sin(y(12))) - (k2./tp).*(y(8)/y(7)).*(sin(y(12))) + (k3./tp).*(y(9)./y(8))*sin(y(13)) + (k1./tp).*(y(6)/y(7)).*sin(y(11));
dy(13) = o3 - (k3./tp).*(y(8)/y(9)).*(sin(y(13))) - (k3./tp).*(y(9)/y(8)).*(sin(y(13))) + (k4./tp).*(y(10)./y(9))*sin(y(14)) + (k2./tp).*(y(7)/y(8)).*sin(y(12));
dy(14) = o4 - (k4./tp).*(y(9)/y(10)).*(sin(y(14))) - (k4./tp).*(y(10)/y(9)).*(sin(y(14))) + (k5./tp).*(y(6)./y(10))*sin(y(15)) + (k3./tp).*(y(8)/y(9)).*sin(y(13));
dy(15) = o5 - (k5./tp).*(y(10)/y(6)).*(sin(y(15))) - (k5./tp).*(y(6)/y(10)).*(sin(y(15))) + (k1./tp).*(y(7)./y(6))*sin(y(11)) + (k4./tp).*(y(9)/y(10)).*sin(y(14));

end

Answer 1

[Sam Chak]

Hi @SAHIL SAHOO

Your k vector is assigned incorrectly. Try if this code works for you.

ti    = 0;
tf    = 1E-3;
tspan = [ti tf];
k     = linspace(1, 2*pi, 5);   % create 5 evenly-spaced points between 1 and 2pi
for i = 1:length(k)
    y0    = [1; 1; 1; 1; 1; 1; 1; 1; 1; 1; k(i); 0; 0; 0; 0];   % just need y0 for each iteration for passing to ode45
    tp    = 5.4E-9;
    [T,Y] = ode45(@(t,y) rate_eq(t, y), tspan, y0);
    E(i)  = mean(Y(:,11));
    q(i)  = (5.*(log(E(i))))./(2*pi);
    plot(T, Y(:,11)), hold on       % just to test if there are multiple plots
end
hold off

function dy = rate_eq(t, y)
    
    dy = zeros(3*5, 1);
    
    % parameters
    P = 1.67;
    a = 0.1;
    tc = 230E-6;
    tp = 5.4E-9;
    o1 = 1E4;
    o2 = 2E4;
    o3 = 3E4;
    o4 = 4E4;
    o5 = 11E4;
    k1 = 5.4E-3;
    k2 = 10.8E-3;
    k3 = 16.2E-3;
    k4 = 22E-3;
    k5 = 54E-3;
    
    % dynamics
    dy(1)  = (P - y(1).*((abs(y(6)))^2 +1))./tc; 
    dy(2)  = (P - y(2).*((abs(y(7)))^2 +1))./tc; 
    dy(3)  = (P - y(3).*((abs(y(8)))^2 +1))./tc; 
    dy(4)  = (P - y(4).*((abs(y(9)))^2 +1))./tc; 
    dy(5)  = (P - y(5).*((abs(y(10)))^2 +1))./tc; 
    dy(6)  = (y(1)-a).*((y(6))./tp) + (k1./tp).*(y(7)).*cos(y(11)) + (k5./tp).*(y(10))*cos(y(15));
    dy(7)  = (y(2)-a).*((y(7))./tp) + (k2./tp).*(y(8)).*cos(y(12)) + (k1./tp).*(y(6))*cos(y(11));
    dy(8)  = (y(3)-a).*((y(8))./tp) + (k3./tp).*(y(9)).*cos(y(13)) + (k2./tp).*(y(7))*cos(y(12));
    dy(9)  = (y(4)-a).*((y(9))./tp) + (k4./tp).*(y(10)).*cos(y(14)) + (k3./tp).*(y(8))*cos(y(13));
    dy(10) = (y(5)-a).*((y(10))./tp) + (k5./tp).*(y(9)).*cos(y(14)) +(k4./tp).* (y(6))*cos(y(15));
    dy(11) = o1 - (k1./tp).*(y(6)/y(7)).*(sin(y(11))) - (k1./tp).*(y(10)/y(6)).*(sin(y(11))) + (k2./tp).*(y(8)./y(7))*sin(y(12)) + (k5./tp).*(y(10)/y(6)).*sin(y(15));
    dy(12) = o2 - (k2./tp).*(y(7)/y(8)).*(sin(y(12))) - (k2./tp).*(y(8)/y(7)).*(sin(y(12))) + (k3./tp).*(y(9)./y(8))*sin(y(13)) + (k1./tp).*(y(6)/y(7)).*sin(y(11));
    dy(13) = o3 - (k3./tp).*(y(8)/y(9)).*(sin(y(13))) - (k3./tp).*(y(9)/y(8)).*(sin(y(13))) + (k4./tp).*(y(10)./y(9))*sin(y(14)) + (k2./tp).*(y(7)/y(8)).*sin(y(12));
    dy(14) = o4 - (k4./tp).*(y(9)/y(10)).*(sin(y(14))) - (k4./tp).*(y(10)/y(9)).*(sin(y(14))) + (k5./tp).*(y(6)./y(10))*sin(y(15)) + (k3./tp).*(y(8)/y(9)).*sin(y(13));
    dy(15) = o5 - (k5./tp).*(y(10)/y(6)).*(sin(y(15))) - (k5./tp).*(y(6)/y(10)).*(sin(y(15))) + (k1./tp).*(y(7)./y(6))*sin(y(11)) + (k4./tp).*(y(9)/y(10)).*sin(y(14));

end