嵌套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双重集成

回答

相关问题