2

我有这段代码：

time = 614.4;
Uhub = 11;
HubHt = 90;
TI = 'A';
N1 = 4096;
N2 = 32;
N3 = 32;

L1 = Uhub*time;
L2 = 150;
L3 = 220;

V = L1*L2*L3;

gamma = 3.9;
c = 1.476;
b = 5.6;

if HubHt < 60
    lambda1 = 0.7*HubHt;
else
    lambda1 = 42;
end
L = 0.8*lambda1;

if isequal(TI,'A')
    Iref = 0.16;
    sigma1 = Iref*(0.75*Uhub + b);
elseif isequal(TI,'B')
    Iref = 0.14;
    sigma1 = Iref*(0.75*Uhub + b);
elseif isequal(TI,'C')
    Iref = 0.12;
    sigma1 = Iref*(0.75*Uhub + b);
else
    sigma1 = str2num(TI)*Uhub/100;
end

sigma_iso = 0.55*sigma1;


%% Wave number vectors
ik1 = cat(2,(-N1/2:-1/2),(1/2:N1/2));
ik2 = -N2/2:N2/2-1;
ik3 = -N3/2:N3/2-1;

[x y z] = ndgrid(ik1,ik2,ik3);

k1 = reshape((2*pi*L/L1)*x,N1*N2*N3,1);
k2 = reshape((2*pi*L/L2)*y,N1*N2*N3,1);
k3 = reshape((2*pi*L/L3)*z,N1*N2*N3,1);

k = sqrt(k1.^2 + k2.^2 + k3.^2);

现在我应该计算

在此处输入图像描述

在哪里

在此处输入图像描述

计算积分的过程是

在此处输入图像描述

目前我正在使用这个循环

E = @(k) (1.453*k.^4)./((1 + k.^2).^(17/6));

E_int = zeros(1,N1*N2*N3);
E_int(1) = 1.5;
for i = 2:(N1*N2*N3) 
    E_int(i) = E_int(i) + quad(E,i-1,i);
end

忽略 k>400 的近似值。我相信我的循环不对。

您建议如何计算积分？

我提前谢谢你。

WKR，弗朗切斯科

4

1 回答 1

2

这是一个从更明显到可能更微妙的更正列表。（事实上​​，我从你在最后一部分写的向上开始）。

从你写的：

 E = @(k) (1.453*k.^4)./((1 + k.^2).^(17/6));

 E_int = zeros(1,N1*N2*N3);
 E_int(1) = 1.5;
 for i = 2:(N1*N2*N3) 
    %//No point in doing this:
    %//E_int(i) = E_int(i) + quad(E,i-1,i); 
    %//According to what you write, it should be:
    E_int(i) = E_int(i-1) + quad(E,i-1,i);     
 end

你可以通过做加速整个事情

%//Independent integration on segments 
Local_int = arrayfun(@(i)quad(E,i-1,i), 2:(N1*N2*N3)); 
Local_int = [1.5  Local_int];
%//integral additivity
E_int = cumsum(Local_int);

此外，如果已知条件（第 2 点）确实是“ ... ( = 1.5 if k' = 0)”，那么整个实现应该更像是

%//Independent integration on segments 
Local_int = arrayfun(@(i)quad(E,i-1,i), 2:(N1*N2*N3));      
%//integral additivity + cumulative removal of queues
E_int = 1.5 - [0 fliplr(cumsum(fliplr(Local_int)))]; %//To remove queues
于 2012-11-24T13:09:21.850 回答