2012-11-06 77 views
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); 

現在我想執行這個積分

enter image description here

與phi_11代替phi_ii。

我該如何在Matlab中執行它,因爲我有所有的矢量化變量?

我提前感謝您的支持。

西九龍填海區, 弗朗切斯科

回答

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); 

您可能需要調整一個整體因爲​​和k2的尺寸會相當隨意(dblquad是一種自適應方法,它將根據函數的本地行爲對點進行評估) 。

+0

您可以在此鏈接中找到hypergeometric2f1函數:http://www.mathworks.com/matlabcentral/fileexchange/35008-generation-of-random-variates/content/pfq.m。不幸的是,目前,我無法訪問Matlab來測試你的建議。 – fpe