Geting error while using quadgk for quadruple integral

I am trying to evaluate quadruple integral using quadgk. It is giving me error even with just one integral. I am pasting the code below

a = 1;
b = 1;
r = 20;
intgrl = quadgk(@(theta1)quadgk(@(theta2)quadgk(@(w)quadgk(@(z)qfunc((a./w.^2)+(b./z.^2)).*(w.^2+z.^2-(4*r^2)<=0).*(w.^2+z.^2-2.*w.*z.*cos(theta1-theta2)>=0).*w.*z.*(1/16*pi^2*r^4), r, 2*r), 0, r), 0,2*pi), 0, 2*pi);

And I am getting the following error.

"Error using + Matrix dimensions must agree."

I feel there is dimension mismatch when evaluating over different variables. Could someone let me know how to evaluate this quadruple integral?

Answer 1

[Mike Hosea]

If you want to integrate a function of 4 variables, f(theta1,theta2,w,z), try iterating it as a single integral of a triple integral like so.

Q1 = @(theta1)integral3(@(theta2,w,z)f(theta1,theta2,w,z),theta2min,theta2max,wmin,wmax,zmin,zmax);

Now you can polish it off by

intgrl = integral(Q1,theta1min,theta1max,'ArrayValued',true);

The 'ArrayValued',true part ensures that integral() will not try to evaluate its integrand function, Q1, at anything other than scalar values of theta1. I'm not saying that this is the fastest way, but the 'ArrayValued' option is the easiest way because otherwise you have to do the work to define your integrand function, f, in such a way that it satisfies the requirements of integral or quadgk, namely that if given arrays of the same size as inputs for all its arguments, it returns an array of the same size. But quadgk returns a scalar, not an array, so you would need to write loops or use arrayfun to define your integrand function, and it can get complicated, or at least complicated-looking.

[Sri]

I have got another problem. When integrating as you suggested, for a particular value of parameter 'r', I keep getting warnings. And for other values of 'r', it is taking very long time to evaluate integral.

Instead, if I use quadv, I am getting results in less time and even for that value of 'r', which gave warnings in the other way. The code is given below

r = 1;
a = 1.13e+04;
b = a;
v = 1;
const = (1/(16*pi^2*r^4));

%using integral
intn = @(theta1) integral3(@(theta2, ds, di) qfunc((a./ds.^v)-(a./di.^v)).*(ds.^2+di.^2-(4*r^2)<=0).*(ds.^2+di.^2-2.*ds.*di.*cos(theta1-theta2)>=0).*ds.*di,r,2*r, 0,r,0,2*pi);

tintn = 0.5*integral(intn, 0,2*pi,'ArrayValued',true)*const

%using quadv
intn = quadv(@(theta1) quadv(@(theta2) quadv(@(ds) quadv(@(di) qfunc((a./ds.^v)-(a./di.^v)).*(ds.^2+di.^2-(4*r^2)<=0).*(ds.^2+di.^2-2.*ds.*di.*cos(theta1-theta2)>=0).*ds.*di,r,2*r), 0,r),0,2*pi),0,2*pi);
tintn = 0.5*intn*const

So, is quadv better to use? Kindly let me know.

[Mike Hosea]

Accuracy-wise, QUADV is not reliable enough for this problem, I think. I don't trust it.

The 'tiled' method of INTEGRAL2 and INTEGRAL3 hates masking functions. In fact masking functions are always a problem numerically, but the 'tiled' method just can't handle them. It has to do with the way it refines in 2D rather than just in 1D. You can switch to the 'iterated' method of INTEGRAL3, though you will probably need to loosen the tolerances substantially for performance reasons. For best results, you should try to reformulate this so that you are integrating qfunc((a./ds.^v)-(a./di.^v))) over the correctly defined domain instead of integrating over a rectangular region and masking it off to be zero over part of the rectangle. The outer boundary of the innermost double integral seems to be the circle ds^2 + di^2 <= (2*r)^2, but the inner boundary is a little more complicated. Bottom line is that defining the region of integration using function handles should work a LOT better if you can sort it all out. It might require splitting the integral up, but it is so expensive to integrate over the discontinuities that the masking functions create, you'll find it to be a bargain.