MATLAB: How to integrate multivariable function numerically

numerical integration for multivariable function with mpower

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?

Best Answer

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
Related Question