matlabでのfor文を減らしたい

Question

下記のプログラムのfor文の使用を1回と0回に減らしたいのですがfor文の使用回数を減らすにはどのように修正すればよいでしょうか？

clear all;
close all;

h=0;
rho=1000;
Ae=0.01;
t=1;
v=1500;
omega=2pi1000000;
k=omega/v;
r=0.1;

x_start=-0.003;
x_end=0.003;
num_div=50;
delta_x=(x_end-x_start)/num_div;
x=x_start:delta_x:x_end;
y_start=-0.003;
y_end=0.003;
delta_y=(y_end-y_start)/num_div;
y=y_start:delta_y:y_end;
[X,Y]=meshgrid(x,y);

P1=0;

for i=0:pi/60:pi/3
for j=0:pi/10:(2*pi-pi/10)

a=i;
b=j;

x1=rcos(b)sin(a);
y1=rsin(b)sin(a);
z1=rcos(a);
deltaS=r^2sin(a)(pi/60)(pi/10);
R=sqrt((X-x1).^2+(Y-y1).^2+(z1).^2);
P=(rhoAe.exp(1iomegat).deltaS)./(2pi.(R+h).exp(1ik(R+h)));

P1=P1+P;

end
end

P1=abs(P1);
mesh(X,Y,P1)

grid on

Comments

[WathMorks]

パラメータ付けされた関数(P)から、パラメータについて和をとることであらたに関数(P1)を定義しているようですが、このような関数をmeshで表示する場合は、for loopを解消しても計算速度は変わらないでしょう。

Answer 1

自己解決しました。3次元配列を使うことで速くできました。
clear all;
close all;
tic
h=0;
rho=1000;
Ae=0.01;
t=1;
v=1500;
omega=2pi1000000;
k=omega/v;

x_start=-0.003;
x_end=0.003;
num_div=50;
delta_x=(x_end-x_start)/num_div;
x=x_start:delta_x:x_end;
y_start=-0.003;
y_end=0.003;
delta_y=(y_end-y_start)/num_div;
y=y_start:delta_y:y_end;

r=0.1;
theta_start=0;
theta_end=pi/3;
num_div=20;
delta_theta=(theta_end-theta_start)/num_div;
theta=theta_start:delta_theta:theta_end;
phi_start=0;
phi_end=2pi-(pi/10);
num_div=19;
delta_phi=(phi_end-phi_start)/num_div;
phi=phi_start:delta_phi:phi_end;
[X,Theta,Phi]=meshgrid(x,theta,phi);
x1=rcos(Phi).sin(Theta);
y1=rsin(Phi).sin(Theta);
z1=rcos(Theta);

P2=zeros(51,51);

for i=1:51

deltaS=r^2sin(Theta)(pi/60)(pi/10);
R=sqrt((X-x1).^2+(y(1,i)-y1).^2+(z1).^2);
P=(rhoAe.exp(1iomegat).deltaS)./(2pi.(R).exp(1ik*(R)));
P1=sum(P);
P1=sum(P1,3);
P2(:,i)=P1;
end

P2=abs(P2);
mesh(P2)
toc

grid on

Answer 2

forを一つにする。（計算速度は2重ループとたいしてかわらないはず）
'''
[a,b]=meshgrid(0:pi/200:pi/3,0:pi/10:(2pi-pi/10));
x1=rcos(b).sin(a);
y1=rsin(b).sin(a);
z1=rcos(a);
deltaS=r^2sin(a)(pi/60)(pi/10);
P1=0;
for j=1:length(a(:))
R=sqrt((X-x1(j)).^2+(Y-y1(j)).^2+(z1(j)).^2);
P=(rhoAeexp(1iomegat)deltaS(j))./(2pi.(R+h).exp(1ik*(R+h)));
P1=P1+P;
end
P1hh=abs(P1);
'''