nested trapz double integration

Question

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

where PHI11 is defined as

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

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

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

Answer 1

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

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.

First, I define test values:

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

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

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

This gives the correct answer.

---

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)

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.

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)

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.