提高5D矩阵计算

Question

早安大家，

我已经写了这个代码：

clc 
clear all 
close all 

%Box size 
Nx=4096; 
Ny=15; 
Nz=15; 

%Spatial gird resolution 
delta=6; 

%WT/turbulence condition 
UHub=11.4; 
HubHt=90; 
z0=0.03; 
IECturbC='B'; 

%%INITIALISATION 

% definition of constants 
twopi=2*pi; 
fourpi=4*pi; 
sqrt2=sqrt(2); 

%constants and derived parameters from IEC 
gamma = 3.9; %IEC, (B.12) 
alpha = 0.2; %IEC, sect. 6.3.1.2 

%set delta1 according to guidelines (chap.6) 
if HubHt<=60, 
    delta1=0.7*HubHt; 
else 
    delta1=42; 
end; 

%IEC, Table 1, p.22 
if IECturbC == 'A', 
    Iref=0.16; 
elseif IECturbC == 'B', 
    Iref=0.14; 
elseif IECturbC == 'C', 
    Iref=0.12; 
else 
    error('IECturbC can be equal to A,B or C;adjust the input value') 
end; 

%IEC, sect. 6.3.1.3 
b=6.5; 
sigma1=Iref*(0.75*UHub+b); 
%derived constants 
l=0.8*delta1; %IEC, (B.12) 
sigmaiso=0.55*sigma1; %IEC, (B.12) 

%%MAIN PROGRAM 
Cij=zeros(3,3,Nx,Ny,Nz); 
k = zeros(3,1); %current vector k 

for ikx=1:(Nx), 
    m = -1.*Nx/2+ikx; 
    k(1)=m*l/(Nx*delta)*twopi; 
    for iky=1:(Ny), 
     m= -1.*Ny/2+iky; 
     k(2)=m*l/(Ny*delta)*twopi; 
     for ikz=1:(Nz), 
     m= -1.*Nz/2+ikz; 
     k(3)=m*l/(Nz*delta)*twopi; 


      if k(1)==0, 
      Cij(:,:,ikx,iky,ikz)=0; 
      else 
      kabs=sqrt(k(1)^2+k(2)^2+k(3)^2); 
      beta= gamma./(kabs.^(2/3)); 
      k0(3)=k(3)+beta.*k(1); 
      k0abs=sqrt(k(1)^2+k(2)^2+k0(3)^2); 
      Ek0=1.453*k0abs^4/(1.+k0abs.^2)^(17/6); 
      C1=beta.*k(1)^2*(k0abs.^2 - 2*k0(3)^2 + beta.*k(1)*k0(3))/(kabs.^2*(k(1)^2 + k(2)^2)); 
      C2=k(2).*k0abs.^2./ (exp((3/2).*log(k(1).^2 + k(2).^2))) .* atan2(beta.*k(1).* sqrt(k(1)^2 + k(2)^2) ,(k0abs.^2 - k0(3).*k(1).*beta)); 
      xhsi1=C1 - k(2).*C2./k(1); 
      xhsi2=k(2).*C1./k(1) + C2; 

      Cij(1,1,ikx,iky,ikz)= sigmaiso*sqrt(twopi*pi*l^3.*Ek0/(Nx*Ny*Nz*delta^3.*k0abs.^4))*(k(2).*xhsi1); 
      Cij(1,2,ikx,iky,ikz)= sigmaiso*sqrt(twopi*pi*l^3.*Ek0/(Nx*Ny*Nz*delta^3.*k0abs.^4))*(k(3) - k(1).*xhsi1 + beta.*k(1)); 
      Cij(1,3,ikx,iky,ikz)= sigmaiso*sqrt(twopi*pi*l^3.*Ek0/(Nx*Ny*Nz*delta^3.*k0abs.^4))*(-k(2)); 
      Cij(2,1,ikx,iky,ikz)= sigmaiso*sqrt(twopi*pi*l^3.*Ek0/(Nx*Ny*Nz*delta^3.*k0abs.^4))*(k(2).*xhsi2 - k(3) - beta.*k(1)); 
      Cij(2,2,ikx,iky,ikz)= sigmaiso*sqrt(twopi*pi*l^3.*Ek0/(Nx*Ny*Nz*delta^3.*k0abs.^4))*(-k(1).*xhsi2); 
      Cij(2,3,ikx,iky,ikz)= sigmaiso*sqrt(twopi*pi*l^3.*Ek0/(Nx*Ny*Nz*delta^3.*k0abs.^4))*(k(1)); 
      Cij(3,1,ikx,iky,ikz)= sigmaiso*sqrt(twopi*pi*l^3.*Ek0/(Nx*Ny*Nz*delta^3.*k0abs.^4))*(k0abs.^2.*k(2) ./ (kabs.^2)); 
      Cij(3,2,ikx,iky,ikz)= sigmaiso*sqrt(twopi*pi*l^3.*Ek0/(Nx*Ny*Nz*delta^3.*k0abs.^4))*(-k0abs.^2*k(1) ./ (kabs.^2)); 
      Cij(3,3,ikx,iky,ikz)= 0; 
      end;  
     end; 
    end; 
end; 

我想问你： 1.有没有更快的方法来获得Cij的矩阵？当Nx，Ny，Nz增加时，Cij的计算相当慢; 2.有什么方法可以获得情节（kabs，beta）和情节（kabs，Ek0）吗？

请耐心等待，我仍然是matlab世界中的新手。

提前感谢和问候， 弗朗切斯科

Answer 1

如果你想在MATLAB良好performaces，你应该尝试矢量化你的代码尽可能。

例如，而不是做：

for x=1:n 
    A(x)=x^2 
end 

做

x=1:n; 
A=x.^2; 

当你有一个以上的指标，您可以使用ndgrid。因此，而不是这样做的：

for x=1:nx 
    for y=1:ny 
    for z=1:nz 
     A(x,y,z)=x^2+y-2*z; 
    end 
    end 
end 

做

[x y z]=ndgrid(1:nx,1:ny,1:nz) 
A=x.^2+y-2*z 

既然你看起来就像你努力，我已经改变了你的代码你。执行时间现在是0.33秒。该vecotrized版本是：

clc 
clear all 
close all 
tic 

%Box size 
Nx=1024; 
Ny=15; 
Nz=15; 

%Spatial gird resolution 
delta=6; 

%WT/turbulence condition 
UHub=11.4; 
HubHt=90; 
z0=0.03; 
IECturbC='B'; 

%%INITIALISATION 

% definition of constants 
twopi=2*pi; 
fourpi=4*pi; 
sqrt2=sqrt(2); 

%constants and derived parameters from IEC 
gamma = 3.9; %IEC, (B.12) 
alpha = 0.2; %IEC, sect. 6.3.1.2 

%set delta1 according to guidelines (chap.6) 
if HubHt<=60, 
    delta1=0.7*HubHt; 
else 
    delta1=42; 
end; 

%IEC, Table 1, p.22 
if IECturbC == 'A', 
    Iref=0.16; 
elseif IECturbC == 'B', 
    Iref=0.14; 
elseif IECturbC == 'C', 
    Iref=0.12; 
else 
    error('IECturbC can be equal to A,B or C;adjust the input value') 
end; 

%IEC, sect. 6.3.1.3 
b=6.5; 
sigma1=Iref*(0.75*UHub+b); 
%derived constants 
l=0.8*delta1; %IEC, (B.12) 
sigmaiso=0.55*sigma1; %IEC, (B.12) 

Cij2=zeros(3,3,Nx,Ny,Nz); 
[x y z]=ndgrid(1:Nx,1:Ny,1:Nz); 
k1=(x-Nx/2)*l/(Nx*delta)*twopi; 
k2=(y-Ny/2)*l/(Ny*delta)*twopi; 
k3=(z-Nz/2)*l/(Nz*delta)*twopi; 
kabs=sqrt(k1.^2+k2.^2+k3.^2); 
beta= gamma./(kabs.^(2/3)); 
k03=k3+beta.*k1; 
k0abs=sqrt(k1.^2+k2.^2+k03.^2); 
Ek0=1.453*k0abs.^4./(1+k0abs.^2).^(17/6); 
C1=beta.*k1.^2.*(k0abs.^2 - 2*k03.^2 + beta.*k1.*k03)./(kabs.^2.*(k1.^2 + k2.^2)); 
C2=k2.*k0abs.^2./ (exp((3/2).*log(k1.^2 + k2.^2))) .* atan2(beta.*k1.* sqrt(k1.^2 + k2.^2) ,(k0abs.^2 - k03.*k1.*beta)); 
xhsi1=C1 - k2.*C2./k1; 
xhsi2=k2.*C1./k1 + C2; 
CC=sigmaiso*sqrt(twopi*pi*l^3.*Ek0./(Nx*Ny*Nz*delta^3.*k0abs.^4)); 
Cij2(1,1,:,:,:)= CC.*(k2.*xhsi1); 
Cij2(1,2,:,:,:)= CC.*(k3 - k1.*xhsi1 + beta.*k1); 
Cij2(1,3,:,:,:)= CC.*(-k2); 
Cij2(2,1,:,:,:)= CC.*(k2.*xhsi2 - k3 - beta.*k1); 
Cij2(2,2,:,:,:)= CC.*(-k1.*xhsi2); 
Cij2(2,3,:,:,:)= CC.*(k1); 
Cij2(3,1,:,:,:)= CC.*(k0abs.^2.*k2 ./ (kabs.^2)); 
Cij2(3,2,:,:,:)= CC.*(-k0abs.^2.*k1 ./ (kabs.^2)); 

Answer 2

我会尽量更普遍地回答你的问题，所以其他人也可以受益。

您应该矢量化for循环以使您的代码更快。取而代之的是这样的：

for i=1:n 
    for j=1:m 
     M(i,j)=sqrt(i) + sqrt(j); 
    end 
end 

矢量化根据下面的代码循环：

[xi,xj] = ndgrid(1:n,1:m); 
M = sqrt(xi)+sqrt(xj);