Using integral3 multiple times

I would like to do numerical triple integration several times and evaluate the final (sextuple) integral at any point I want. For example, let's say I'd like to do the following:

a = -2.5;
b = 2.5;
f = @(x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8) exp(-(x_1.*x_2.*x_3.*x_4.*x_5.*x_6.*x_7.*x_8));
g = @(x_1,x_2,x_3,x_4,x_5) integral3(@(x_6,x_7,x_8) f(x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8) , a,b, a,b,a,b);
k = @(x_1,x_2) integral3(@(x_3,x_4,x_5) g(x_1,x_2,x_3,x_4,x_5) , a,b, a,b,a,b);
k_num=k(0,0)

then I get the following message:

Error using .* Matrix dimensions must agree.

Do you have any idea how to fix it?

Answer 1

[Steven Lord]

From the documentation:

"The function fun must accept three arrays of the same size and return an array of corresponding values. It must perform element-wise operations."

There's no guarantee that the integral3 call inside the anonymous function k passes the same size arrays as x_3, x_4, and x_5 into g (and because g calls f, into f) that the integral3 call inside g passes as x_6, x_7, and x_8 into f.

But as written your problem is easy. You call k with x_1 equal to 0. k calls g with x_1 equal to 0. g calls f with x_1 equal to 0. f multiplies x_1 through x_8 together (resulting in 0), multiplies that by -1 (again resulting in 0) and then computes exp(0) which equals 1.

Since I'm guessing this was just for demonstration purposes and your real problem is more complicated, g will need to integrate f once per element in the x_3, x_4, and x_5 arrays that were passed into it by k. You could use arrayfun in the anonymous function to do that. This isn't going to be fast because it's going to call integral3 many times, but performing a sextuple integral isn't that common an operation in my experience so there is no integral6.

[Apostolos Psaros]

Thanks! You are right, my function is more complicated. I would like to try arrayfun though, just to see how fast it is.

I'm doing the following but it doesn't work:

f = @(x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8) exp(-(x_1.*x_2.*x_3.*x_4.*x_5.*x_6.*x_7.*x_8));
fv = @(x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8) arrayfun(f,x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8);

a = -2.5;
b = 2.5;

g = @(x_1,x_2,x_3,x_4,x_5) integral3(@(x_6,x_7,x_8) fv(x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8) , a,b, a,b,a,b);

g(0,0,0,0,0)

Could you see the mistake, if any?