decreasing three loops to two loops with sum

1 visualización (últimos 30 días)
Roya Rajabi
Roya Rajabi el 23 de Sept. de 2020
Comentada: Roya Rajabi el 23 de Sept. de 2020
I have this code and want to eliminate the inner most loop and use sum
How can I do that?
thetao=1;
alfa=1;
l=1;
theta2(1:7,1:21)=0;
T=[0.01, 0.1, 0.2, 0.5, 1, 2, 5];
x=0:0.05:l;
n=0:50;
for t=1:7
for j=1:length(x)
for k=1:length(n)
theta1(t,j)=(sin((2*n(k)+1)*pi*x(j)/(2*l))*exp(-(2*n(k)+1).^2*pi^2*alfa*T(t)/(4*l^2))/((2*n(k)+1)*pi));
theta2(t,j)=theta2(t,j)+theta1(t,j);
theta(t,j)=thetao-4*thetao*theta2(t,j);
end
end
end

Respuesta aceptada

Adam Danz
Adam Danz el 23 de Sept. de 2020
Editada: Adam Danz el 23 de Sept. de 2020
----------- Updated answer ------------------------------
Vectorize the k-loop and use sum() to define theta2.
thetao=1;
alfa=1;
l=1;
theta2(1:7,1:21)=0;
T=[0.01, 0.1, 0.2, 0.5, 1, 2, 5];
x=0:0.05:l;
n=0:50;
for t=1:7
for j=1:length(x)
theta1 = sin((2*n+1)*pi*x(j)/(2*l)).* exp(-(2*n+1).^ 2*pi^2*alfa*T(t)/(4*l^2))./((2*n+1)*pi);
theta2 = sum(theta1);
theta(t,j)=thetao-4*thetao*theta2;
end
end
To demonstrate that this produces the same values at the original version, run the block and above and then two two blocks below.
% Save the theta values above and then clear all variables
theta_version2 = theta;
clearvars -except theta_version2
% Original version
thetao=1;
alfa=1;
l=1;
theta2(1:7,1:21)=0;
T=[0.01, 0.1, 0.2, 0.5, 1, 2, 5];
x=0:0.05:l;
n=0:50;
for t=1:7
for j=1:length(x)
for k=1:length(n)
theta1(t,j)=(sin((2*n(k)+1)*pi*x(j)/(2*l))*exp(-(2*n(k)+1).^2*pi^2*alfa*T(t)/(4*l^2))/((2*n(k)+1)*pi));
theta2(t,j)=theta2(t,j)+theta1(t,j);
theta(t,j)=thetao-4*thetao*theta2(t,j);
end
end
end
Now compare the results
max(abs(theta - theta_version2),[],'all')
% ans =
% 4.4409e-16
The max difference in theta between the original and new versions is 4.4409e-16. This is due to the precision in floating point calculations explained [here].
  3 comentarios
Adam Danz
Adam Danz el 23 de Sept. de 2020
See my updated answer.
Roya Rajabi
Roya Rajabi el 23 de Sept. de 2020
Thanks a million it works

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