is there any method to apply the for loop in this, I caN't see any pattern please help in this

2 views (last 30 days)
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

Accepted Answer

David Goodmanson
David Goodmanson on 28 Aug 2022
Edited: David Goodmanson on 28 Aug 2022
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 Comments
Torsten
Torsten on 29 Aug 2022
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));

Sign in to comment.

More Answers (0)

Categories

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