How to integrate multivariable function numerically?

12 次查看(过去 30 天)
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?

采纳的回答

John D'Errico
John D'Errico 2016-1-12
编辑:John D'Errico 2016-1-12
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

更多回答(1 个)

jgg
jgg 2016-1-12
编辑:jgg 2016-1-12
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

类别

Help CenterFile Exchange 中查找有关 Linear Algebra 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by