How to integrate multivariable function numerically?

I need to perform numerical integration for

exp(x^2+y^2+z^2)

When I use

 f = exp(x^2+y^2+z^2);
 F = @(x,y,z)(exp(x^2+y^2+z^2));
 Q = triplequad(F,0,1,0,1,0,1)

It gives ??? Error using ==> mpower Matrix must be square. Is there a way to perform this?

Answer 1

[John D'Errico]

Of course, this being a separable problem, it is trivially doable, in far less time than triplequad will take. As jgg showed, the problem that had to be resolved to use triplequad was to make your function properly vectorized. Thus, use .^ instead of ^ as an operator. Now it can use inputs of any shape or size, as long as x, y, and z are all the same shape and size. You would also need to use .* and ./ if multiplication and division were used.

Next, MATLAB tells us that triplequad is being replaced with integral3.

format long g
tic
F = @(x,y,z)(exp(x.^2+y.^2+z.^2));
Q = integral3(F,0,1,0,1,0,1)
toc
Q =
        3.12912420249805

Elapsed time is 0.065742 seconds.

tic
F = @(x,y,z)(exp(x.^2+y.^2+z.^2));
Q = triplequad(F,0,1,0,1,0,1)
toc
Q =
        3.12912428413994

Elapsed time is 0.172366 seconds.

So integral 3 is a bit faster than is triplequad anyway. As it turns out, integral3 was more accurate too.

The point is though, since it is separable, just write the integrand as:

   exp(x^2)*exp(y^2)*exp(z^2)

and then recognize that the domain of integration is a simple rectangular one. So the problem really is separable.

format long g
tic,I1 = integral(f,0,1);toc
Elapsed time is 0.002503 seconds.

I1^3
ans =
     3.12912420246652

I'll take the latter, especially since it will be more accurate yet.

Is that claim true? Lets see, since we can do a bit better though yet, since each individual integral in the separable problem is solvable analytically.

For general upper and lower limits, we get a result in the form of the special function erfi.

syms x a b
Ix = int(exp(x^2),a,b)
Ix =
-(pi^(1/2)*(erfi(a) - erfi(b)))/2

And with limits of [0,1], we get:

subs(Ix,{a,b},{0,1})
ans =
(pi^(1/2)*erfi(1))/2

finally, since the limits are the same for each integral, we get the result

format long g
((pi^(1 / 2)*erfi(1))/2)^3
ans =
     3.12912420246652

It looks like the separable result was indeed quite accurate, since we can compare these results to a symbolic one:

vpa('((pi^(1 / 2)*erfi(1))/2)^3')
ans =
3.1291242024665164843069648833865

Answer 2

[jgg]

The issue is that for triplequad to work fun(x,y,z) must accept a vector x and scalars y and z, and return a vector of values of the integrand. Your function doesn't do this. I think you just need this:

 F = @(x,y,z)(exp(x.^2+y.^2+z.^2));

So that it computes a vector of inputs. Then, proceed as before:

 Q = triplequad(F,0,1,0,1,0,1) %3.1291 as expected