is there any method to apply the for loop in this, I caN't see any pattern please help in this
1 次查看(过去 30 天)
显示 更早的评论
ti = 0;
tf = 10E-2;
tspan=[ti tf];
KC = 1E-3;
y0= ones(1,80)*1E-3;
yita_mn = [
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0;
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1;
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
]*(KC);
N = 20;
tp = 5.4E-9;
[T,Y]= ode45(@(t,y) rate_eq(t,y,yita_mn,N),tspan,y0);
% i want to use for loop in this
plot(T./tp,(Y(:,61)));
hold on
plot(T./tp,(Y(:,61)));
plot(T./tp,(Y(:,62)));
plot(T./tp,(Y(:,63)));
plot(T./tp,(Y(:,64)));
plot(T./tp,(Y(:,65)));
plot(T./tp,(Y(:,66)));
plot(T./tp,(Y(:,67)));
plot(T./tp,(Y(:,68)));
plot(T./tp,(Y(:,69)));
plot(T./tp,(Y(:,70)));
plot(T./tp,(Y(:,71)));
plot(T./tp,(Y(:,72)));
plot(T./tp,(Y(:,73)));
plot(T./tp,(Y(:,74)));
plot(T./tp,(Y(:,75)));
plot(T./tp,(Y(:,76)));
plot(T./tp,(Y(:,77)));
plot(T./tp,(Y(:,78)));
plot(T./tp,(Y(:,79)));
plot(T./tp,(Y(:,80)));
hold off
function dy = rate_eq(t,y,yita_mn,N)
dy = zeros(4*N,1);
dGdt = zeros(N,1);
dAdt = zeros(N,1);
dOdt = zeros(N,1);
P = 16.7;
a = 0.1;
tc = 230E-6;
tp = 5.4E-9;
o = rand(1,20).*10E3;
Gt = y(1:3:3*N-2);
At = y(2:3:3*N-1);
Ot = y(3:3:3*N-0);
k = 0.0001;
for i = 1:N
dGdt(i) = (P - Gt(i).*((abs(At(i)))^2 +1))./tc ;
dAdt(i) = (Gt(i)-a).*((At(i))./tp);
dOdt(i) = o(1,i);
for j = 1:N
dAdt(i) = dAdt(i)+yita_mn(i,j)*(1/tp)*(At(j))*cos(Ot(j)-Ot(i));
dOdt(i) = dOdt(i)+yita_mn(i,j)*(1/tp)*((At(j)/At(i)))*sin(Ot(j)-Ot(i));
end
end
dy(1:3:3*N-2) = dGdt;
dy(2:3:3*N-1) = dAdt;
dy(3:3:3*N-0) = dOdt;
%here i want to compress this. but couldn't understand how to apply for loop here.
dy(61) = o(1,2) - o(1,1) - (k./tp).*(y(2)/y(5)).*(sin(y(61))) - (k./tp).*(y(5)/y(2)).*(sin(y(61))) + (k./tp).*(y(8)./y(5))*sin(y(62)) + (k./tp).*(y(59)/y(2)).*sin(y(80));
dy(62) = o(1,3) - o(1,2) - (k./tp).*(y(5)/y(8)).*(sin(y(62))) - (k./tp).*(y(8)/y(5)).*(sin(y(62))) + (k./tp).*(y(11)./y(8))*sin(y(63)) + (k./tp).*(y(2)/y(5)).*sin(y(61));
dy(63) = o(1,4) - o(1,3) - (k./tp).*(y(8)/y(11)).*(sin(y(63))) - (k./tp).*(y(11)/y(8)).*(sin(y(63))) + (k./tp).*(y(14)./y(11))*sin(y(64)) + (k./tp).*(y(5)/y(8)).*sin(y(62));
dy(64) = o(1,5) - o(1,4) - (k./tp).*(y(11)/y(14)).*(sin(y(64))) - (k./tp).*(y(14)/y(11)).*(sin(y(64))) + (k./tp).*(y(17)./y(14))*sin(y(65)) + (k./tp).*(y(8)/y(11)).*sin(y(63));
dy(65) = o(1,6) - o(1,5) - (k./tp).*(y(14)/y(17)).*(sin(y(65))) - (k./tp).*(y(17)/y(14)).*(sin(y(65))) + (k./tp).*(y(20)./y(17))*sin(y(66)) + (k./tp).*(y(11)/y(14)).*sin(y(64));
dy(66) = o(1,7) - o(1,6) - (k./tp).*(y(17)/y(20)).*(sin(y(66))) - (k./tp).*(y(20)/y(17)).*(sin(y(66))) + (k./tp).*(y(23)./y(20))*sin(y(67)) + (k./tp).*(y(14)/y(17)).*sin(y(65));
dy(67) = o(1,8) - o(1,7) - (k./tp).*(y(20)/y(23)).*(sin(y(67))) - (k./tp).*(y(23)/y(20)).*(sin(y(67))) + (k./tp).*(y(26)./y(23))*sin(y(68)) + (k./tp).*(y(17)/y(20)).*sin(y(66));
dy(68) = o(1,9) - o(1,8) - (k./tp).*(y(23)/y(26)).*(sin(y(68))) - (k./tp).*(y(26)/y(23)).*(sin(y(68))) + (k./tp).*(y(29)./y(26))*sin(y(69)) + (k./tp).*(y(20)/y(23)).*sin(y(67));
dy(69) = o(1,10) - o(1,9) - (k./tp).*(y(26)/y(29)).*(sin(y(69))) - (k./tp).*(y(29)/y(26)).*(sin(y(69))) + (k./tp).*(y(32)./y(29))*sin(y(70)) + (k./tp).*(y(23)/y(26)).*sin(y(68));
dy(70) = o(1,11) - o(1,10) - (k./tp).*(y(29)/y(32)).*(sin(y(70))) - (k./tp).*(y(32)/y(29)).*(sin(y(70))) + (k./tp).*(y(35)./y(32))*sin(y(71)) + (k./tp).*(y(26)/y(29)).*sin(y(69));
dy(71) = o(1,12) - o(1,11) - (k./tp).*(y(32)/y(35)).*(sin(y(71))) - (k./tp).*(y(35)/y(32)).*(sin(y(71))) + (k./tp).*(y(38)./y(35))*sin(y(72)) + (k./tp).*(y(29)/y(32)).*sin(y(70));
dy(72) = o(1,13) - o(1,12) - (k./tp).*(y(35)/y(38)).*(sin(y(72))) - (k./tp).*(y(38)/y(35)).*(sin(y(72))) + (k./tp).*(y(41)./y(38))*sin(y(73)) + (k./tp).*(y(32)/y(35)).*sin(y(71));
dy(73) = o(1,14) - o(1,13) - (k./tp).*(y(38)/y(41)).*(sin(y(73))) - (k./tp).*(y(41)/y(38)).*(sin(y(73))) + (k./tp).*(y(44)./y(41))*sin(y(74)) + (k./tp).*(y(35)/y(38)).*sin(y(72));
dy(74) = o(1,15) - o(1,14) - (k./tp).*(y(41)/y(44)).*(sin(y(74))) - (k./tp).*(y(44)/y(41)).*(sin(y(74))) + (k./tp).*(y(47)./y(44))*sin(y(75)) + (k./tp).*(y(38)/y(41)).*sin(y(73));
dy(75) = o(1,16) - o(1,15) - (k./tp).*(y(44)/y(47)).*(sin(y(75))) - (k./tp).*(y(47)/y(44)).*(sin(y(75))) + (k./tp).*(y(50)./y(47))*sin(y(76)) + (k./tp).*(y(41)/y(44)).*sin(y(74));
dy(76) = o(1,17) - o(1,16) - (k./tp).*(y(47)/y(50)).*(sin(y(76))) - (k./tp).*(y(50)/y(47)).*(sin(y(76))) + (k./tp).*(y(53)./y(50))*sin(y(77)) + (k./tp).*(y(44)/y(47)).*sin(y(75));
dy(77) = o(1,18) - o(1,17) - (k./tp).*(y(50)/y(53)).*(sin(y(77))) - (k./tp).*(y(53)/y(50)).*(sin(y(77))) + (k./tp).*(y(56)./y(53))*sin(y(78)) + (k./tp).*(y(47)/y(50)).*sin(y(76));
dy(78) = o(1,19) - o(1,18) - (k./tp).*(y(53)/y(56)).*(sin(y(78))) - (k./tp).*(y(56)/y(53)).*(sin(y(78))) + (k./tp).*(y(59)./y(56))*sin(y(79)) + (k./tp).*(y(50)/y(53)).*sin(y(77));
dy(79) = o(1,20) - o(1,19) - (k./tp).*(y(56)/y(59)).*(sin(y(79))) - (k./tp).*(y(59)/y(56)).*(sin(y(79))) + (k./tp).*(y(2)./y(59))*sin(y(80)) + (k./tp).*(y(53)/y(56)).*sin(y(78));
dy(80) = o(1,1) - o(1,20) - (k./tp).*(y(57)/y(2)).*(sin(y(80))) - (k./tp).*(y(2)/y(57)).*(sin(y(80))) + (k./tp).*(y(5)./y(2))*sin(y(61)) + (k./tp).*(y(56)/y(59)).*sin(y(79));
end
0 个评论
采纳的回答
David Goodmanson
2022-8-28
编辑:David Goodmanson
2022-8-28
Hi Sahil,
This can be done by for loop but it is maybe better to do the whole thing with index manipulation. In the written-out expressions for dy, looking down the columns of indices there are some circular shifts between some of the columns. There are two kinds of indices, those that increment by 1, and the y indices that increment by 3. Denoting the first set of indices with n, then
n1 = (1:20)';
n2 = circshift(n1,-1);
n61 = n1 +60;
n62 = circshift(n61,-1);
n80 = circshift(n61,1);
% as a check, compare the indices for o( ) and sin(y( )) with the written-out version
nset = [n1 n2 n61 n62 n80]
nset =
1 2 61 62 80
2 3 62 63 61
3 4 63 64 62
4 5 64 65 63
5 6 65 66 64
6 7 66 67 65
7 8 67 68 66
8 9 68 69 67
9 10 69 70 68
10 11 70 71 69
11 12 71 72 70
12 13 72 73 71
13 14 73 74 72
14 15 74 75 73
15 16 75 76 74
16 17 76 77 75
17 18 77 78 76
18 19 78 79 77
19 20 79 80 78
20 1 80 61 79
Here n80 is called that because its value is 80 in the first row, and similarly for the others, first row.
For the increment-by-3 indices,
j2 = 3*(1:20)-1;
j5 = circshift(j2,-1);
j8 = circshift(j2,-2);
j59 = circshift(j2,1);
% as a check, compare these y( ) indices with the written-out version
jset = [j2; j5; j8; j59]'
jset =
2 5 8 59
5 8 11 2
8 11 14 5
11 14 17 8
14 17 20 11
17 20 23 14
20 23 26 17
23 26 29 20
26 29 32 23
29 32 35 26
32 35 38 29
35 38 41 32
38 41 44 35
41 44 47 38
44 47 50 41
47 50 53 44
50 53 56 47
53 56 59 50
56 59 2 53
59 2 5 56
NOTE that your expression for dy(80) contains two values of 57 which by all rights look like they should be 59.
Given the indices, all of which are vectors of length 20, you can calculate all of the dy in one go:
dy(n61) = o(1,n2) - o(1,n1) ...
- (k./tp).*(y(j2) ./y(j5)).*sin(y(n61)) - (k./tp).*(y( j5) ./y(j2)).*sin(y(n61)) ...
+ (k./tp).*(y(j8) ./y(j5)).*sin(y(n62)) + (k./tp).*(y(j59) ./y(j2)).*sin(y(n80));
for the plotting:
figure(1)
plot(T./tp,(Y(:,61)));
hold on
for m = 62:80
plot(T./tp,(Y(:,m)));
end
hold off
2 个评论
Torsten
2022-8-29
dy(n61) = o(1,n2).' - o(1,n1).' - (k./tp).*(y(j2)./y(j5)).*sin(y(n61)) - (k./tp).*(y( j5)./y(j2)).*sin(y(n61)) + (k./tp).*(y(j8)./y(j5)).*sin(y(n62)) + (k./tp).*(y(j59)./y(j2)).*sin(y(n80));
instead of
dy(n61) = o(1,n2) - o(1,n1) - (k./tp).*(y(j2)./y(j5)).*sin(y(n61)) - (k./tp).*(y( j5)./y(j2)).*sin(y(n61)) + (k./tp).*(y(j8)./y(j5)).*sin(y(n62)) + (k./tp).*(y(j59)./y(j2)).*sin(y(n80));
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!