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

1 visualización (últimos 30 días)
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
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
SAHIL SAHOO
SAHIL SAHOO el 29 de Ag. de 2022
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);
Unable to perform assignment because the left and right sides have a different number of elements.

Error in solution>rate_eq (line 80)
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));

Error in solution (line 30)
[T,Y]= ode45(@(t,y) rate_eq(t,y,yita_mn,N),tspan,y0);

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
figure(1)
plot(T./tp,(Y(:,61)));
hold on
for m = 62:80
plot(T./tp,(Y(:,m)));
end
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 = [6 2.3 4.4 5.3 4.6 7.8 5.4 3.5 4.6 7.3 8.7 9.5 5.6 6.6 7.6 4.5 3.4 2.5 3 5.3 4.5]*10E3;
Gt = y(1:3:3*N-2);
At = y(2:3:3*N-1);
Ot = y(3:3:3*N-0);
k = 0.001;
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;
n1 = (1:20)';
n2 = circshift(n1,-1);
n61 = n1 +60;
n62 = circshift(n61,-1);
n80 = circshift(n61,1);
j2 = 3*(1:20)-1;
j5 = circshift(j2,-1);
j8 = circshift(j2,-2);
j59 = circshift(j2,1);
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));
end
i apply your comments in this still this is not working, what the problem in this.
Torsten
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));

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Loops and Conditional Statements en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by