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_11代替phi_ii。
我該如何在Matlab中執行它,因爲我有所有的矢量化變量?
我提前感謝您的支持。
西九龍填海區, 弗朗切斯科
您可以在此鏈接中找到hypergeometric2f1函數:http://www.mathworks.com/matlabcentral/fileexchange/35008-generation-of-random-variates/content/pfq.m。不幸的是,目前,我無法訪問Matlab來測試你的建議。 – fpe