嵌套trapz双重集成 [英] nested trapz double integration

问题描述

我想知道是否可以通过嵌套的trapz循环绕过对quad2d的调用.我将更详细地讨论我的问题:目前，我以这种方式执行双积分的计算:

I'd like knowing if there is any way to bypass the call for quad2d with a nested trapz loop. I'll discuss my problem with some more detail: currently, I perform the calculation of a double integral this way:

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可以通过以下方式获得:

and E_integral.mat can be obtained this way:

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

现在的问题是:通过使用嵌套的trapz函数，是否有可能忽略quad2d和handle function而采用更实用的方法?

Now the question is: would it be possible to overlook quad2d and the handle function in favor of a more practicle approach, by using a nested trapz function?

到目前为止，我已经尝试了以下代码，但未获得预期的结果:

So far, I've tried the following piece of code, which has not yielded to the expected results:

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

推荐答案

我不确定为什么要避免使用quad2d函数，但是可以将2D正交分解为两个1D嵌套的正交(即使该函数不会分解).

I'm not sure why you want to avoid using the quad2d function, but it's possible to split the 2D quadrature into two 1D nested quadratures (even if the function does not factorize).

这是通过两个trapz调用进行集成的一种处理方式.关键是要在一个大表中生成所有值的集合，并在两个维度上使用trapz.

Here is a way to process with the integration with two trapz calls. The point is to generate the set of all the values in a big table and the use trapz for both dimensions.

首先，我定义测试值:

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

然后，我同时定义x和y方向的网格(类似于您的k2和k3):

Then, I define the grid for both the x and y directions (similar to your k2 and k3):

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

然后我评估所有点对中的函数f:

I evaluate then the function f in all the pairs of points:

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

然后，可以使用两个嵌套调用来完成集成:

Then, the integration can be done using two nested calls:

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

这给出了正确的答案.

在我的答案的第一个版本中，我使用了quad函数，该函数直接与函数句柄一起使用. 天真"的想法是简单地嵌套两个调用，例如:

In the first version of my answer, I had used the quad function, which works with function handles directly. The "naive" idea is to simply nest the two calls, like:

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

，但这实际上不起作用.关键是内部的quad调用不适用于矢量值的函​​数，而外部的调用则需要矢量值的函​​数.

but this does not work actually. The point is that the inner quad call does not work for vector valued functions while the outer call expects a vector values function.

要使其正常工作，我们只需要更改quad，quadv的向量值版本的内部调用:

To make it work, we simply need to change the inner call for the vector valued version of quad, quadv:

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

现在它应该可以工作了，我们可以检查一下它给出的答案是否正确

Now it should work and we can check that it actually gives the same answer as

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

如果您真的想使用trapz函数，建议您在trapz之上构建一个函数，该函数接受函数句柄并在内部调用trapz.

If you really want to use the trapz function, I suggest that you build a function on top of trapz that accepts a function handle and calls trapz inside.

这篇关于嵌套trapz双重集成的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！