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??

1 view (last 30 days)
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

Accepted Answer

Sam Chak
Sam Chak on 22 Aug 2022
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

More Answers (0)

Categories

Find more on Historical Contests 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!