执行分步积分

Question

我有一段代码：

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，弗朗西斯科

Answer 1

这是一个从更明显的到可能更微妙的更正列表。(实际上，我从您在最后部分中写的内容开始往上看)。

从你写的来看：

 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