在Matlab中使用向量化循环进行双重求和
|
我想对这个double for循环进行矢量化处理，因为这是我代码中的瓶颈。由于Matlab是一种基于一种的索引语言，因此我必须为M = 0创建一个附加项。
R，r，lambda是常数
Slm（L，M），（L，M）是矩阵70x70
Plm（L，M）是70x71的矩阵
Cl（L），Pl（L）是向量70x1
% function dU_r
s1 = 0;
for L = 2:70
s1 = s1 + ((R/r)^L)*(L+1)*Pl(L)*Cl(L);
for m = 1:L
s1 = s1 + ((R/r)^L)*(L+1)*Plm(L,M)*(Clm(L,M)*...
cos(M*lambda) + Slm(L,M)*sin(M*lambda));
end
end
dU_r = -(mu_p/(r^2))*s1;
% function dU_phi
s2=0;
for L = 2:70
s2 = s2 + ((R/r)^L))*Plm(L,1)*Cl(L);
for m = 1:l
s2 = s2 + ((R/r)^L)*(Plm(L,M+1)-M*tan(phi)*Plm(L,M))*...
(Clm(L,M)*cos(M*lambda) + Slm(L,M)*sin(M*lambda));
end;
end;
dU_phi = (mu_p/r)*s2;
% function dU_lambda
s3=0;
for L=2:70
for m=1:L
s3 = s3 + ((R/r)^L)*M*Plm(L,M)*(Slm(L,M)*cos(M*lambda)...
- Clm(L,M)*sin(M*lambda));
end;
dU_lambda = (mu_p/r)*s3;
没有找到相关结果
已邀请:
3 个回复
坛沤疲撑拆
，，和/或。s2 = 0; L = 2:70; PlCl = Pl.*Cl; PlmClm = Plm.*Clm; Rrl = (R/r).^(L).*(L+1); COS = cos((1:70)*lambda); SIN = sin((1:70)*lambda); for l=L M = 1:l; s1 = sum(Rrl(l-1)*PlmClm(l,M).*(COS(M) + Slm(l,M).*SIN(M))); s2 = s2 + s1; end; s2 = s2 + sum(Rrl(L).*PlCl(L));辅奈
，和，则仅使用它们较低三角形部分的值。主对角线上的所有值都将被忽略。通过对向量进行完全矢量化，使其使用矩阵运算而不是循环，您最终将对矩阵的上三角部分进行清零，从而执行不必要的乘法。这是for循环可能是更好选择的情况之一。 因此，这是Egon答案的精炼和更正版本，就执行的数学运算次数而言，应该接近最佳值：index = 2:70; coeff = ((R/r).^index).*(index+1); lambda = lambda.*(1:70).\'; %\' cosVector = cos(lambda); sinVector = sin(lambda); s2 = coeff*(Pl(index).*Cl(index)); for L = index M = 1:L; s2 = s2 + coeff(L-1)*((Plm(L,M).*Clm(L,M))*cosVector(M) + ... (Plm(L,M).*Slm(L,M))*sinVector(M)); end炬卤遁蝎变
Rr = R/r; RrL = Rr; cosLambda = cos((1:70)* lambda); sinLambda = sin((1:70)* lambda); u1 = uint8(1); s1 = 0; s2 = 0; s3 = 0; for j = uint8(2):uint8(70) RrL = RrL * Rr; q1 = RrL * (double(j) + 1); t1 = Pl(j) * datos.Cl(j); q2 = RrL; t2 = Plm(j,1) * Cl(j); t3 = 0; for m = u1:j t1 = t1 + Plm(j,m) * ... (Clm(j, m) * cosLambda(m) + ... Slm(j, m) * sinLambda(m)); t2 = t2 + (Plm(j,m+1)-double(m)*tan_phi*Plm(j,m))*... (Clm(j,m)*cosLambda(m) + Slm(j,m)*sinLambda(m)); t3 = t3 + double(m)*Plm(j,m)*(Slm(j,m)*cosLambda(m)... - Clm(j,m)*sinLambda(m)); end s1 = s1 + q1 * t1; s2 = s2 + q2 * t2; s3 = s3 + q2 * t3; end dU_r = -(mu_p/(r^2))*s1; dU_phi = (mu_p/r)*s2; dU_lambda = (mu_p/r)*s3