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

Question

SAHIL SAHOO 2022-8-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1783165-in-this-i-apply-a-loop-so-that-i-will-have-different-value-of-y-11-for-the-different-value-of-ini

回答： Sam Chak 2022-8-22

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Sam Chak 2022-8-22

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1783165-in-this-i-apply-a-loop-so-that-i-will-have-different-value-of-y-11-for-the-different-value-of-ini#answer_1030475

在 MATLAB Online 中打开

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