2 views (last 30 days)
Sri on 8 Jul 2014
Commented: Mike Hosea on 2 Sep 2014
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;
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?

Mike Hosea on 11 Aug 2014
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.
##### 2 CommentsShowHide 1 older comment
Mike Hosea on 2 Sep 2014
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.