1

答えるのは簡単な質問があると思います。

私はこのコードを持っています：

clc, clear all, close all

Nx = 8192;
Ny = 16;
Nz = 16;
delta = 6;
UHub = 11.4;
IECturbC = 'A';
HubHt = 90;

%%INITIALISATION

% definition of constants
twopi=2*pi;
%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=5.6;
sigma1=Iref*(0.75*UHub+b);
sigma2=0.7*sigma1;
sigma3=0.5*sigma1;
%derived constants
l=0.8*delta1; %IEC, (B.12)
sigmaiso=0.55*sigma1; %IEC, (B.12)

Cij2=zeros(3,3,Nx,Ny,Nz);
ikx=cat(2,-Nx/2:-1,1:Nx/2-1);
[x y z]=ndgrid(ikx,1:Ny,1:Nz);
k1=(x)*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).*(hypergeometric2f1(1/3, 17/6,4/3,-kabs.^(-2),11))).^(0.5));
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));

次に、私は定義します

phi_11 = (Ek0./4*pi.*kabs.^4).*(k0abs.^2 - k1.^2 - 2*k1.*k03.*xhsi1 + (k1.^2 + k2.^2).*xhsi1.^2);

今、私はこの積分を実行したいと思います

ここに画像の説明を入力してください

phi_iiの代わりにphi_11を使用します。

すべてのベクトル化された変数があるので、Matlabでどのように実行できますか？

よろしくお願いします。

WKR、フランチェスコ

4

1 に答える 1

2

私はhypergeometric2f1（関数？変数？）を持っていないので、これをテストすることはできません。

とにかく、これを実行する通常の方法は、次のphi_11ように定義することです。

phi_11 = @(k1,k2) (Ek0./4*pi.*kabs.^4).*(k0abs.^2 - k1.^2 - 2*k1.*k03.*xhsi1 + (k1.^2 +  k2.^2).*xhsi1.^2);

次に、二重求積法を使用します。

F_11 = dblquad(phi_11, -inf,inf, -inf,inf);

k1とはかなり任意の次元にk2なるため、変数全体を調整する必要がある場合があります（dblquadこれは適応型の方法であり、関数のローカル動作に依存するポイントで関数を評価します）。

于 2012-11-06T09:01:24.510 に答える