numerical integration - Quadruple Integral Using Nested Integral2 in Matlab

Question

I am trying to solve a problem of the following form:

f=@(x,y,z,w) x.*y.*z.*w;          % A complicated black box function
a=1;b=1;c=1;d=1;                  % Integration limits
I=integral2(@(x,y)integral2(@(z,w)f(x,y,z,w),c,-c,d,-d),a,-a,b,-b);

Using this implementation I get the following error:

Error using  .* 
Matrix dimensions must agree.

The problem is that x, y, z, and w are not the same size. For the first function evaluation all inputs are the same size but then on the second function evaluation x and y are not the same size as z and w.

How can I resolve this error?

This question is similar to this unanswered question: Input array size error for a quadraple integration using nested integral2

==================================================================================

In response to the answer:

I=integral(@(x)integral3(@(y,z,w)f(x,y,z,w),b,-b,c,-c,d,-d),a,-a,'ArrayValued',true);

This does solve the problem, however, it is not obvious to me why this works. I had actually seen this solution before but forgot to mention it in my question (http://www.mathworks.com/matlabcentral/answers/77571-how-to-perform-4d-integral-in-matlab).

I would like to solve using nested integral2 because I know my function is discontinuous and would like to use the iterated integration method. I could do something like this but only the inner integral is iterated so I am not sure how that affects accuracy:

I=integral(@(x)integral3(@(y,z,w)f(x,y,z,w),b,-b,c,-c,d,-d,'Method','iterated'),a,-a,'ArrayValued',true);

Answer 1

integral2 calls its integrand with two matrix arguments of the same size. The problem is that you can't just mix the variables x,y,z,w in the function call f(x,y,z,w) as you did in your question, because the dimension of x and y is determined by the outer integral2, whereas the dimension of z and w is determined by the inner integral2, so is isn't guaranteed that the dimensions are the same. The integralX functions don't vectorize anyway, each call can provide just one output value.

The function integral provides an option so that it calls its integrand only with scalar values, and with Matlab scalar expansion, this works together with the same-sized 3D arrays the inner function integral3 provides.

I=integral(@(x)integral3(@(y,z,w)f(x,y,z,w),b,-b,c,-c,d,-d),a,-a,'ArrayValued',true);

You can achieve the same (calling with scalars) by encapsulating the inner integral2 call with arrayfun:

I=integral2(@(x,y)arrayfun(@(x,y)integral2(@(z,w)f(x,y,z,w),c1,c2,d1,d2),x,y),a1,a2,b1,b2)

The latter was about six times faster in my experiments.