how to calculate the Guass function integral within a ellipsoid?

function:F=exp(-(2*(x^2+y^2))/a^2-2*z^2/b^2)

integral within a ellipsoid space:(x^2+y^2)/a^2+z^2/b^2<1 a>0,b>0 Thanks a lot!

Answer 1

[Teja Muppirala]

The INTEGRAL3 function can do this quite easily. For example, for a = 2 and b = 3,

a = 2;
b = 3;

f = @(x,y,z) exp(-(2*(x.^2+y.^2))/a^2-2*z.^2/b^2);
xmin = -a;
xmax = a;
ymin = @(x)-sqrt(a^2 - x.^2);
ymax = @(x) sqrt(a^2 - x.^2);
zmin = @(x,y)-sqrt(b^2*(1-(x.^2+y.^2)/a^2));
zmax = @(x,y) sqrt(b^2*(1-(x.^2+y.^2)/a^2));
Q = integral3(f,xmin,xmax,ymin,ymax,zmin,zmax)

You get about 17.44 as the answer.

Just as a sanity check, you could also do it the old fashioned way by making a grid and summing things up:

N = 200;
[x,y,z] = ndgrid(linspace(-a,a,N),linspace(-a,a,N),linspace(-b,b,N));
dx = x(2,1,1) - x(1,1,1);
dy = y(1,2,1) - y(1,1,1);
dz = z(1,1,2) - z(1,1,1);

CONSTRAINT = ( (x.^2+y.^2)/a^2+z.^2/b^2 < 1 );
F = exp(-(2*(x.^2+y.^2))/a^2-2*z.^2/b^2) .* CONSTRAINT;
Q = sum(F(:))*dx*dy*dz

Indeed, this gives you pretty much the same answer as INTEGRAL3 does.

[Henry Liu]

Thanks very much for your illuminating answer！

However，I want to get the results of symbolic integration. Can you help me？

[Teja Muppirala]

Oh. Well symbolically, you should rescale things and do it in spherical coordinates, and it's pretty simple.

Define x' = x/a, y' = y/a z' = z/b

Then your integral

int [ exp(-(2*(x^2+y^2))/a^2-2*z^2/b^2) dx dy dz]

becomes

a^2*b*int [ exp(-(2*(x'^2+y'^2+z'^2))) dx' dy' dz']

You have basically scaled the ellipsoid so that you are integrating over a unit sphere instead.

Then you can change to spherical coordinates:

a^2*b*int [ exp(-(2*(r^2))) r^2 sin(theta) dr dtheta dphi]

Where r = 0..1, theta = 0..pi phi = 0..2*pi

Since theta and phi are not in the integrand, you can bring them out in front as a 4*pi,

Then we just need to calculate

4*pi*a^2*b*int [ exp(-(2*(r^2))) r^2 dr] for r = 0..1

Putting it all together, in MATLAB you get:

syms a b r
Q = simplify( 4*pi*a^2*b*int(exp(-(2*(r^2)))*r^2,0,1) )
double( subs(Q,{a,b},{2,3}) )

You can see that the analytic expression obtained is in agreement with the numeric answer from before.

[Henry Liu]

Wonderful！ Thanks very much！