decreasing three loops to two loops with sum

1 回表示 (過去 30 日間)
Roya Rajabi
Roya Rajabi 2020 年 9 月 23 日
コメント済み: Roya Rajabi 2020 年 9 月 23 日
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

採用された回答

Adam Danz
Adam Danz 2020 年 9 月 23 日
編集済み: Adam Danz 2020 年 9 月 23 日
----------- 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 件のコメント
Adam Danz
Adam Danz 2020 年 9 月 23 日
See my updated answer.
Roya Rajabi
Roya Rajabi 2020 年 9 月 23 日
Thanks a million it works

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeLoops and Conditional Statements についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by