Multiple integrals

Definition

Multiple integrals are definite integrals and they arise in many areas of physics, in particular, in mechanics, where volumes, masses, and moments of inertia of bodies are of interest. Definition of multiple integrals is the extention of the definition of simple integrals: The volume of integration V is being split into many small cubes or other small elements of the volumes (Subscript[ΔV, i]) and the integral is defined as the limit of the sum. For a 3d integral one defines

(∫_V)f[r]dV≡∫∫(∫_V)f[x,y,z]dxdydz=(lim_(N→∞))(∑_(i=1)^N)f[(Subscript[r, i])](Subscript[ΔV, i])

In particular, for the moment of inertia around the z axis (Subscript[I, zz]) one has

(Subscript[I, zz]) =ρ(∫_V)((x^2)+(y^2))dV≡ρ∫∫(∫_V)((x^2)+(y^2))dxdydz=(lim_(N→∞))(∑_(i=1)^N)(Subscript[r, i,⊥]^2)(Subscript[ΔM, i])

where (Subscript[ΔM, i]) = ρ (Subscript[ΔV, i]) and ρ is the mass density.

Recursive calculation of integrals

Another equivalent definition of a multiple integral is recursive and reduces it to simple integrals. In particular, a double integral can be introduced as the repeated integral of the type

∫∫f[x,y]dxdy=(∫_(Subscript[x, 1])^(Subscript[x, 2]))F[x]dx,      F[x]=(∫_((Subscript[y, 1])[x])^((Subscript[y, 2])[x]))f[x,y]dy.

It can be shown that the result does not depend on the order of integrations over x and y and is the same as according to the symmetric definition above. As an example les us calculate the area of a rectangle -a/2 < x < a/2,  -b/2 < y < b/2:

S=∫(∫_(-a/2 < x < a/2
-b/2 < y < b/2))dxdy=(∫_(-a/2)^(a/2))dx (∫_(-b/2)^(b/2))dy=(∫_(-a/2)^(a/2))dx b = ab

The same can be done by the Mathematica command

(∫_(-b/2)^(b/2))(∫_(-a/2)^(a/2))\[DifferentialD]x\[DifferentialD]y

or

(∫_(-a/2)^(a/2))(∫_(-b/2)^(b/2))\[DifferentialD]y\[DifferentialD]x

or

Integrate[1,{x,-a/2,a/2},{y,-b/2,b/2}]

Another example is the area of a circle r < R.

(Subscript[S=∫(∫_(r < R))dxdy=(∫_(-R)^R)dx (∫_(-(Sqrt[(R^2)-(x^2)]))^(Sqrt[(R^2)-(x^2)]))dy=(∫_(-R)^R)dx 2(Sqrt[(R^2)-(x^2)])=(x (Sqrt[(R^2)-(x^2)])+(R^2) ArcTan[(x/(Sqrt[(R^2)-(x^2)]))])|, -R]^R)=((πR)^2)

In this case, there is no elegant Mathematica command to integrate over the region (x^2)+(y^2)<(R^2). Instead, one has to do work and set the integration limits explicitly

(∫_(-R)^R)(∫_(-(Sqrt[(R^2)-(x^2)]))^(Sqrt[(R^2)-(x^2)]))\[DifferentialD]y\[DifferentialD]x //Timing

or

Integrate[1,{x,-R,R},{y,-(Sqrt[(R^2)-(x^2)]),(Sqrt[(R^2)-(x^2)])},Assumptions->R>0] //Timing

Ellipsoid is defined by the equation

((x^2)/(a^2))+((y^2)/(b^2))+((z^2)/(c^2))=1

The volume of the ellipsoid is given by the integral

Vellipsoid=Integrate[1,{x,-a,a},{y,-b(Sqrt[1-((x^2)/(a^2))]),b(Sqrt[1-((x^2)/(a^2))])},{z,-c(Sqrt[1-((x^2)/(a^2))-((y^2)/(b^2))]),c(Sqrt[1-((x^2)/(a^2))-((y^2)/(b^2))])},Assumptions->{a>0,b>0,c>0}]

For the sphere one obtains

Vsphere=Vellipsoid/.{a->R,b->R,c->R}

Moments of inertia of the cone

Let us calculate the momentz of inertia of a cone with the base of radius R and height h. Putting the origin of the coordinate system in the apex of the cone, the z axis being directed from the apex to the base of the cone, one obtains

Izz=ρ Integrate[(x^2)+(y^2),{z,0,h},{x,-(z/h)R,(z/h)R},{y,-(Sqrt[((((z/h)R))^2)-(x^2)]),(Sqrt[((((z/h)R))^2)-(x^2)])},Assumptions->{h>0,R>0}]

Ixx=ρ Integrate[(y^2)+(z^2),{z,0,h},{x,-(z/h)R,(z/h)R},{y,-(Sqrt[((((z/h)R))^2)-(x^2)]),(Sqrt[((((z/h)R))^2)-(x^2)])},Assumptions->{h>0,R>0}]

(Subscript[I, yy])=(Subscript[I, xx]) by symmetry. It is convenient to express the density ρ via the mass of the cone M using ρ = M/V. With

V=Integrate[1,{z,0,h},{x,-(z/h)R,(z/h)R},{y,-(Sqrt[((((z/h)R))^2)-(x^2)]),(Sqrt[((((z/h)R))^2)-(x^2)])},Assumptions->{h>0,R>0}]

one obtains, finally

Izz=Izz/.ρ->(M/V)
Ixx=Ixx/.ρ->(M/V)

The moments of inertia of the cone above are defined with respect to its apex. Also of interest are moments of inertia with respect to the center of mass (Subscript[r, CM]) defined by

(Subscript[r, CM])=(1/V)(∫_V)rdV,     ρ=const

The center of mass of the cone is located, by symmetry, on its z axis,

zCM=(1/V)Integrate[z,{z,0,h},{x,-(z/h)R,(z/h)R},{y,-(Sqrt[((((z/h)R))^2)-(x^2)]),(Sqrt[((((z/h)R))^2)-(x^2)])},Assumptions->{h>0,R>0}]

Now putting the origin of the coordinate system in the center of mass, one obtains another value of the xx moment of inertia, IxxCM

IxxCM=(M/V) Integrate[(y^2)+(((z-zCM))^2),{z,0,h},{x,-(z/h)R,(z/h)R},{y,-(Sqrt[((((z/h)R))^2)-(x^2)]),(Sqrt[((((z/h)R))^2)-(x^2)])},Assumptions->{h>0,R>0}]

Let us check the relation

(Subscript[I, xx])=((Subscript[I, xx])^((CM)))+(Ma^2),

where a is the distance between the center of mass and the point used in the definition of (Subscript[I, xx])

a=zCM;
Ixx==IxxCM +M (a^2) //Simplify

Numerical calculation of multiple integrals

Numerical calculation of multiple integrals in Mathematica is done by the similar command as for their analytical calculation. For instance, consider the trips integral P[G] that appears in the theory of magnetism

λ[x_,y_,z_]=(1/3)(Cos[x]+Cos[y]+Cos[z]);
P[G_]:=(1/(((2π))^3))(∫_(-π)^π)(∫_(-π)^π)(∫_(-π)^π)(1/(1-G λ[x,y,z]))\[DifferentialD]x\[DifferentialD]y\[DifferentialD]z

This and related integrals can be calculated analytically but Mathematica cannot do it because analytical calculation requires a dosen of nontrivial transformations. The command

N[P[1]]

Does not produce any result because Mathematica attempts first to calculate the integral analytically and then find its numerical value at G = 1. To proceed, one has to define this function as a numeric integral

P[G_]:=(1/(((2π))^3))NIntegrate[(1/(1-G λ[x,y,z])),{x,-π,π},{y,-π,π},{z,-π,π}]

P[1]  is the famous Watson integral W. In this case the integrand becomes divergent at x = y = z = 0

Normal[Series[1-λ[x,y,z],{x,0,2},{y,0,2},{z,0,2}]]

This is why Mathematica complains numerically calculating this integral and the accuracy is lower than usual, although not bad.

P[1]

Consulting the Help and increasing the order of the Mutidimensional rule solves the accuracy problem

P[G_]:=(1/(((2π))^3))NIntegrate[(1/(1-G λ[x,y,z])),{x,-π,π},{y,-π,π},{z,-π,π},Method->{"MultiDimensionalRule","Generators"->6}]
 P[1]

Also in this case one can analytically integrate over one of the variables and then integrate over the two remaining variables without problems and faster

PP[G_,x_,y_]=(1/(2π))Integrate[(1/(1-G λ[x,y,z])),{z,-π,π},GenerateConditions->False]

P[G_]:=(1/(((2π))^2))NIntegrate[PP[G,x,y],{x,-π,π},{y,-π,π}]
P[1]

Change of variables in multiple integrals

Polar, spherical, or cylindrical coordinates

If the integration region has a circular, spherical, or cylindrical symmetry, it is convenient to use polar, spherical, or cylindrical coordinates.

Polar coordinates

In two dimensions, one can use the polar coordinates (r,ϕ), instead of the Descarde cordinates (x,y). The relation betwen the two reads

x=r Cos[ϕ]
y=r Sin[ϕ]

The infinitesimal area element dS in the polar coordinates is a nearly-square with sides dr and rdϕ, that is

dS=rdrdϕ

Now two-dimensional integrals can be rewritten as

∫(∫_S)f[x,y]dxdy=∫(∫_S)f[r Cos[ϕ],r Sin[ϕ]]rdrdϕ,

with the appropriate change of integration limits. In particular, the area of a circle is given by a simple integral

Scircle=(∫_0^R)(∫_0^(2π))r\[DifferentialD]ϕ\[DifferentialD]r

or, in the factorized form,

Scircle=((∫_0^R)r\[DifferentialD]r)((∫_0^(2π))\[DifferentialD]ϕ)

The moments of inertia of a circle around the z and x axes are

(Subscript[I, zz])=ρ∫(∫_S)((x^2)+(y^2))dxdy=ρ∫(∫_S)(r^3)drdϕ
(Subscript[I, xx])=ρ∫(∫_S)(y^2)dxdy=ρ∫(∫_S)(r^3)((Sin[ϕ])^2)drdϕ

where ρ is the surface mass density. Integration and substitution ρ = M/S yields

Izz=(M/Scircle)((∫_0^R)(r^3)\[DifferentialD]r)2π
Ixx=(M/Scircle)((∫_0^R)(r^3)\[DifferentialD]r)((∫_0^(2π))((Sin[ϕ])^2)\[DifferentialD]ϕ)

As an application of the polar coordinates, let us consider the trick of calculating the Gaussian integral

I=(∫_(-∞)^∞)(e^(-(x^2)))dx.

Square of this integral can be represented as the product of two similar integrals over x and y, a double integral that can be rewriten in the polar coordinates and easily calculated

(Subscript[(I^2)=((∫_(-∞)^∞)(e^(-(x^2)))dx)((∫_(-∞)^∞)(e^(-(y^2)))dy)=(∫_(-∞)^∞)(∫_(-∞)^∞)(e^(-((x^2)+(y^2))))dxdy=(∫_0^∞)(∫_0^(2π))(e^(-(r^2)))rdϕdr=2π(∫_0^∞)(e^(-(r^2)))rdr=π(∫_0^∞)(e^(-(r^2)))(dr^2)=-((πe)^(-(r^2)))|, 0]^∞)=π.

Thus

I=(Sqrt[π]).

Cylindrical coordinates

Cylindrical coordinates is the straightforward extention of the polar coordinats into 3 dimensions leaving the z axis untouched.

x=r Cos[ϕ]
y=r Sin[ϕ]
z=z

The infinitesimal volume in cylindrical coordinates reads

dV=rdrdϕdz

Now three-dimensional integrals can be rewritten as

∫∫(∫_V)f[x,y,z]dxdydz=∫∫(∫_V)f[r Cos[ϕ],r Sin[ϕ],z]rdrdϕdz

Moment if inertia (Subscript[I, zz]) of a cylinder is given by the same formula as that of a disc. Let us calculate (Subscript[I, xx]) with respect to the center of mass of a cylinder with radius R and height h.

Vcylinder=π (R^2)h;
Ixx=(M/Vcylinder)(∫_(-h/2)^(h/2))(∫_0^R)(∫_0^(2π))((r^2) ((Sin[ϕ])^2)+(z^2))r\[DifferentialD]ϕ\[DifferentialD]r\[DifferentialD]z

For a very short cilinder, h≪R, the result for the disc is reproduced,

Ixx/.h->0

For a very long cilinder, h≫R, the result for the thin rod is reproduced,

Ixx/.R->0

Cylindrical coordinate system is convenient to calculate the volume of a thorus because the thorus has a rotational symmetry around the z axis and integrarion over ϕ yields simply 2π.

R=3;   r=1;
x[u_,v_]:=(R+r Cos[v])Sin[u];
y[u_,v_]:=(R+r Cos[v])Cos[u];
z[u_,v_]:=r Sin[v];
ParametricPlot3D[{x[u,v],y[u,v],z[u,v]},{u,0,2Pi},{v,0,2Pi},PlotStyle->{Red}]

The remaining integration is over a circle of radius r centered at the distance R from the z axis.

R=3;   r=1; z=0.5;
Show[Plot[Null,{ρ,0,R+2r},PlotRange->{-1.5,1.5},AxesLabel->{"ρ","z"}],
Graphics[{Yellow,Rectangle[{0,-r},{R,r}]}],
Graphics[{Orange,Disk[{R,0},r]}],
Graphics[Arrow[{{R-(Sqrt[(r^2)-(z^2)]),z},{R+(Sqrt[(r^2)-(z^2)]),z}}]],
Graphics[Arrow[{{R,-r},{R,r}}]],
Graphics[Disk[{R,z},0.03]]
]

Vthorus=2π(∫_(-r)^r)(∫_(R-(Sqrt[(r^2)-(z^2)]))^(R+(Sqrt[(r^2)-(z^2)])))ρ\[DifferentialD]ρ\[DifferentialD]z

or better

Vthorus=2π Integrate[ρ,{z,-r,r},{ρ,R-(Sqrt[(r^2)-(z^2)]),R+(Sqrt[(r^2)-(z^2)])},Assumptions->{R>0,r>0}]

It is interesting that one can straighten the thorus replacing it by the cylinder of radius r and height 2πR. The volume remains the same and given by V = (((πr)^2))(2πR) = 2(π^2)(r^2) R. The moment of inertia of the thorus with respect to the z axis is given by

Izz=(M/Vthorus)2π Integrate[(ρ^3),{z,-r,r},{ρ,R-(Sqrt[(r^2)-(z^2)]),R+(Sqrt[(r^2)-(z^2)])},Assumptions->{R>0,r>0}]//Simplify

In the limit r = 0 this reduces to the moment of inertia of a ring (MR^2).

Spherical coordinates

Spherical coordinates (r,ϑ,ϕ) are defined by

x=r Sin[ϑ]Cos[ϕ]
y=r Sin[ϑ]Sin[ϕ]
z=r Cos[ϑ]

where ϑ is the polar angle and ϕ is the azimuthal angle. The infinitesimal volume in the spherical coordinates is a nearly-cube with the sides dr, rdϑ, r Sin[ϑ]dϕ, that is,

dV=(r^2)Sin[ϑ]drdϑdϕ

As an illustration, let us calculate the volume of a sphere

Vsphere=(∫_0^R)(∫_0^π)(∫_0^(2π))(r^2)Sin[ϑ]\[DifferentialD]ϕ\[DifferentialD]ϑ\[DifferentialD]r

or, in the factorized form

Vsphere=((∫_0^R)(r^2)\[DifferentialD]r)((∫_0^π)Sin[ϑ]\[DifferentialD]ϑ)((∫_0^(2π))\[DifferentialD]ϕ)

This integration is much simpler than that for the ellipsoid in the Descarde coordinates above. On the oher hand, integration in the spherical coordinates is simple only for a sphere, but for an ellipsoid it becomes complicated.

Let us calculate the moment of inertia of a solid sphere in the spherical coordinates. To find (Subscript[I, zz]), one has to integrate (x^2) + (y^2) = (r^2)((Sin[ϑ])^2) over the sphere's volume

Izz=(M/Vsphere)((∫_0^R)(r^4)\[DifferentialD]r)((∫_0^π)((Sin[ϑ])^3)\[DifferentialD]ϑ)((∫_0^(2π))\[DifferentialD]ϕ)

General change of variables

General change of variables (x,y,z)->(u,v,w) etc in multi-dimensional integrals can be effectuated by

∫∫∫f[x,y,z]dxdydz=∫∫∫f[x[u,v,w],y[u,v,w],z[u,v,w]]((∂(x,y,z))/(∂(u,v,w)))dudvdw,

where J[u, v, w]≡((∂(x,y,z))/(∂(u,v,w))) is the Jacobian of the transformation defined as a determinant of partial derivatives

((∂(x,y,z))/(∂(u,v,w)))≡|({{((∂x)/(∂u)), ((∂x)/(∂v)), ((∂x)/(∂w))}, {((∂y)/(∂u)), ((∂y)/(∂v)), ((∂y)/(∂w))}, {((∂z)/(∂u)), ((∂z)/(∂v)), ((∂z)/(∂w))}})|

In fact the Jacobian is determined as the absolute value of this. The results above are valid in all dimentions. For instance, the Jacobian of the transition to the polar coordinate system is

((∂(x,y))/(∂(r,ϕ)))≡|({{((∂x)/(∂r)), ((∂x)/(∂ϕ))}, {((∂y)/(∂r)), ((∂y)/(∂ϕ))}})|=|({{Cos[ϕ], -r Sin[ϕ]}, {Sin[ϕ], r Cos[ϕ]}})|=r( ((Cos[ϕ])^2)+((Sin[ϕ])^2))=r

in accordance with our previous consideration. Let us calculate the Jacobian of the spherical coordinate system with the help of Mathematica.

x=r Sin[ϑ]Cos[ϕ];
y=r Sin[ϑ]Sin[ϕ];
z=r Cos[ϑ];

J=Det[(({{(∂_r)x, (∂_ϑ)x, (∂_ϕ)x}, {(∂_r)y, (∂_ϑ)y, (∂_ϕ)y}, {(∂_r)z, (∂_ϑ)z, (∂_ϕ)z}}))]//Simplify

again in accordance with previous considerations.

As an illustration of a general change of variables, let us calculate the double integral

II=(∫_0^1)(∫_0^1)(1/(1-(x^2)(y^2)))\[DifferentialD]x\[DifferentialD]y

In terms of new variables (u,v) so that

x=((Sin[u])/(Cos[v]));     y=((Sin[v])/(Cos[u]));

the denominator of the integrand becomes

1-((Tan[u])^2)((Tan[v])^2)

On the other hand, the Jacobian is

J=Det[(({{(∂_u)x, (∂_v)x}, {(∂_u)y, (∂_v)y}}))]//Simplify

Thus in the new variables the integrand is 1. It remains to figure out the integration region in (u,v). Let us first invert the transformation

Solve[{x==((Sin[u])/(Cos[v])), y==((Sin[v])/(Cos[u]))},{u,v}]

Take the positive solution

vxy[x_,y_]:=ArcSec[((Sqrt[-1+(x^2) (y^2)])/(Sqrt[-1+(y^2)]))]
uxy[x_,y_]:=ArcSin[((x (Sqrt[-1+(y^2)]))/(Sqrt[-1+(x^2) (y^2)]))]

Draw the four lines in the (u,v) plane that corresponds to the four sides of the square in the (x,y) plane. As our transformation becomes singular at x = 1 and y = 1, go a bit away from these lines.

a=0.00000000001;
Show[
ParametricPlot[{uxy[x,0],vxy[x,0]},{x,0,1-a},PlotRange->All,PlotStyle->{Thick,Green}],
ParametricPlot[{uxy[x,1-a],vxy[x,1-a]},{x,0,1-a},PlotRange->All,PlotStyle->{Thick,Red}],
ParametricPlot[{uxy[0,y],vxy[0,y]},{y,0,1-a},PlotRange->All,PlotStyle->{Thick,Black}],
ParametricPlot[{uxy[1-a,y],vxy[1-a,y]},{y,0,1-a},PlotRange->All,PlotStyle->{Thick,Blue}]
]

Thus the integration region is bordered by the lines

u=0,    v=0,    u+v=(π/2).

The integral is equal to the area of this triangle,

II=(1/2)×(π/2)×(π/2)=((π^2)/8)

Let us now express our integral via series expanding the integrand

II=(∫_0^1)(∫_0^1)(1/(1-(x^2)(y^2)))dxdy=(∫_0^1)(∫_0^1)(1+(x^2)(y^2)+((((x^2)(y^2)))^2)+((((x^2)(y^2)))^3)+...)dxdy=1+(1/(3^2))+(1/(5^2))+(1/(7^2))+...

This allows to calculate the series

S=(∑_(n=1)^∞)(1/(n^2))

The sum of odd terms of this series is (Subscript[S, odd])=(( (π^2))/8). The sum of its even terms is

(Subscript[S, even])=(1/(2^2))+(1/(4^2))+(1/(6^2))+...=(1/4)(1+(1/(2^2))+(1/(3^2))+...)=(1/4)S

Thus

S=(Subscript[S, even])+(Subscript[S, odd])=(1/4)S+(( (π^2))/8)

and  S=((π^2)/6).