is there any method to apply the for loop in this, I caN't see any pattern please help in this
Mostrar comentarios más antiguos
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
Respuesta aceptada
David Goodmanson
el 28 de Ag. de 2022
Editada: David Goodmanson
el 28 de Ag. de 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 comentarios
Torsten
el 29 de Ag. de 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));
Más respuestas (0)
Categorías
Más información sobre Numerical Integration and Differential Equations en Centro de ayuda y File Exchange.
el 28 de Ag. de 2022
el 29 de Ag. de 2022
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
