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 に答える