嵌套trapz雙重集成

我想知道是否有任何方法繞過呼叫與嵌套trapz循環。我會用一些更詳細的討論我的問題：目前，我執行二重積分這樣計算：嵌套trapz雙重集成

clc, clear all, close all 
load E_integral.mat 

c = 1.476; 
gamma = 3.0; 

beta_int = (c*gamma)./(k_int.*sqrt(E_integral)); 

figure, loglog(k_int,beta_int,'r','LineWidth',2.0), grid on; 

k1 = (.01:.1:100); 
k2 = .01:.1:100; 
k3 = -100:.1:100; 

int_k3 = zeros(size(k2)); 
int_k3k2 = zeros(size(k1)); 

tic 
for ii = 1:numel(k1) 
    phi11 = @(k2,k3) PHI11(k1(ii),k2,k3,k_int,beta_int); 
    int_k3(ii) = 2*quad2d(phi11,-100,100,-100,100); 
end 
toc

其中PHI11被定義爲

function phi11 = PHI11(k1,k2,k3,k_int,beta_int) 
k = sqrt(k1.^2 + k2.^2 + k3.^2); 
ksq = k.^2; 
k1sq = k1.^2; 
fourpi = 4.*pi; 
beta = exp(interp1(log(k_int),log(beta_int),log(k),'linear')); 
k30 = k3 + beta.*k1; 
k0 = sqrt(k1.^2 + k2.^2 + k30.^2); 
k0sq = k0.^2; 
k04sq = k0.^4; 
Ek0 = (1.453.*k04sq)./((1 + k0sq).^(17/6)); 

C1 = (beta.*k1sq.*(k0sq - 2.*(k30.^2) + beta.*k1.*k30))./(ksq.*(k1.^2 + k2.^2)); 
C2 = ((k2.*k0sq)./((k1.^2 + k2.^2).^(3/2))).*atan2((beta.*k1.*sqrt(k1.^2 + k2.^2)),(k0sq - k30.*k1.*beta)); 
xhsi1 = C1 - (k2./k1).*C2; 
xhsi1_sq = xhsi1.^2; 
phi11 = (Ek0./(fourpi.*k04sq)).*(k0sq - k1sq - 2.*k1.*k30.*xhsi1 + (k1.^2 + k2.^2).*xhsi1_sq); 
end

和E_integral.mat可以通過這種方式獲得：

clc,clear all,close all 

k_int = .001:.01:1000; 

Ek = (1.453.*k_int.^4)./((1 + k_int.^2).^(17/6)); 


E = @(k_int) (1.453.*k_int.^4)./((1 + k_int.^2).^(17/6)); 

E_integral = zeros(size(k_int)); 

for ii = 1:numel(k_int)  
    E_integral(ii) = integral(E,k_int(ii),Inf);  
end 

save('E_integral','k_int','E_integral')

現在的問題是：有可能忽略quad2d和最愛的handle function或更實際的方法，通過使用嵌套的trapz函數？

到目前爲止，我已經試過下面的代碼段，這還沒有產生預期的結果：

int_k33 = zeros(size(k2)); 
S_11 = zeros(size(k1)); 
tic 
for ii = 1:1 
    for jj = 1:numel(k2) 
     int_k33(jj) = trapz(k3,PHI11(k1(ii),k2(jj),k3,k_int,beta_int));   
    end 
    S_11(ii) = 4*trapz(k2,int_k33); 
end 
toc

來源

2013-05-06 fpe

任何有趣的想法？ – fpe 2013-05-07 16:34:52

請解釋一下「實用方法」是什麼意思？我只給出了這種前瞻性的外觀，我當然可以看到改進代碼的地方。這就是說，一般來說，如果你可以寫'f（x），那麼你只能將一個雙積分∫∫f（x，y）dA'分解成兩個獨立的積分∫g（x）dx *∫h（y）dy' （x，y）'作爲兩個函數g（x）* h（y）'的乘積。你能發表實際的方程式和他們的上下文嗎？（因子1 /（4pi）使我想到與磁場有關的事情。） – 2013-05-13 21:45:47

我*假設*你想這樣做，因爲你給了實驗數據？否則，選擇這樣一種簡單的方法而不是更精確的方法是沒有意義的......你能詳細說明*爲什麼*你想這樣做？ – 2013-05-16 11:35:52

我不知道爲什麼你要避免使用quad2d功能，但它的可能將二維正交分解爲兩個1D嵌套正交（即使函數沒有因式分解）。

這是一種處理方法，與兩個trapz調用進行集成。重點是生成一個大表中的所有值的集合，並使用兩個維度的使用trapz。

首先，我定義的測試值：

f= @(x,y) sin(x.*cos(y)); 
N = 1000;

然後，我定義爲x和y方向上的網格（類似於您k2和k3）：

xpts1d = linspace(0,1,N+1); 
ypts1d = linspace(0,1,N+1);

然後我評估所有對點中的函數f：

xpts = xpts1d'*ones(1,N+1); 
ypts = ones(N+1,1)*ypts1d; 
values = f(xpts,ypts);

然後，整合日糧可以使用兩個嵌套調用完成：

I = trapz(xpts1d,trapz(ypts1d,values,2),1);

這給出了正確的答案。

在我回答的第一個版本，我已經使用了quad函數，該函數的功能直接處理工作。「天真」的想法是簡單地嵌套兩個電話，如：

I = quad(@(y) quad(@(x)f(x,y) ,0,1),0,1)

但這實際上不起作用。重點在於內部的quad調用不適用於向量值函數，而外部調用需要向量值函數。

要使其工作，我們只需要改變內喊出了quad的向量值版本，quadv：

I = quad(@(y) quadv(@(x)f(x,y) ,0,1),0,1)

現在它應該工作，我們可以檢查它實際上給出了相同的答案

I = quad2d(f, 0,1,0,1)

如果你真的想使用trapz功能，我建議你建立一個接受功能手柄，並呼籲trapz INSI上trapz頂部的功能德。

來源

2013-05-14 13:45:58

嵌套trapz雙重集成

回答

相關問題