Inmathematics (specificallymultivariable calculus), amultiple integral is adefinite integral of afunction of several real variables, for instance,f(x,y) orf(x,y,z).

Integrals of a function of two variables over a region in${\displaystyle {\mathbb{R}}^2}$ (thereal-number plane) are calleddouble integrals, and integrals of a function of three variables over a region in${\displaystyle {\mathbb{R}}^3}$ (real-number 3D space) are calledtriple integrals.^[1] For repeated antidifferentiation of a single-variable function, see theCauchy formula for repeated integration.

Just as the definite integral of a positive function of one variable represents thearea of the region between the graph of the function and thex-axis, thedouble integral of a positive function of two variables represents thevolume of the region between the surface defined by the function (on the three-dimensionalCartesian plane wherez =f(x,y)) and the plane which contains itsdomain.^[1] If there are more variables, a multiple integral will yieldhypervolumes of multidimensional functions.

Multiple integration of a function inn variables:f(x₁,x₂, ...,x_n) over a domainD is most commonly represented by nested integral signs in the reverse order of execution (the leftmost integral sign is computed last), followed by the function and integrand arguments in proper order (the integral with respect to the rightmost argument is computed last). The domain of integration is either represented symbolically for every argument over each integral sign, or is abbreviated by a variable at the rightmost integral sign:^[2]

${\displaystyle \int \cdots {\int }_{\mathbf{D}}f(x_1,x_2,\dots ,x_n)d{x}_1\cdots d{x}_n}$ 

Since the concept of anantiderivative is only defined for functions of a single real variable, the usual definition of theindefinite integral does not immediately extend to the multiple integral.

Forn > 1, consider a so-called "half-open"n-dimensionalhyperrectangular domainT, defined as

${\displaystyle T=[a_1,b_1)\times [a_2,b_2)\times \cdots \times [a_n,b_n)\subseteq {\mathbb{R}}^n}$ .

Partition each interval[a_j,b_j) into a finite familyI_j of non-overlapping subintervalsi_jα, with each subinterval closed at the left end, and open at the right end.

Then the finite family of subrectanglesC given by

${\displaystyle C=I_1\times I_2\times \cdots \times I_n}$ 

is apartition ofT; that is, the subrectanglesC_k are non-overlapping and their union isT.

Letf :T →R be a function defined onT. Consider a partitionC ofT as defined above, such thatC is a family ofm subrectanglesC_m and

${\displaystyle T=C_1\cup C_2\cup \cdots \cup C_m}$ 

We can approximate the total(n + 1)-dimensional volume bounded below by then-dimensional hyperrectangleT and above by then-dimensional graph off with the followingRiemann sum:

${\displaystyle \sum _{k=1}^{m}f(P_k)\mathrm{m}(C_k)}$ 

whereP_k is a point inC_k andm(C_k) is the product of the lengths of the intervals whose Cartesian product isC_k, also known as the measure ofC_k.

Thediameter of a subrectangleC_k is the largest of the lengths of the intervals whoseCartesian product isC_k. The diameter of a given partition ofT is defined as the largest of the diameters of the subrectangles in the partition. Intuitively, as the diameter of the partitionC is restricted smaller and smaller, the number of subrectanglesm gets larger, and the measurem(C_k) of each subrectangle grows smaller. The functionf is said to beRiemann integrable if thelimit

${\displaystyle S=\underset{\delta \to 0}{lim}\sum _{k=1}^{m}f(P_k)\mathrm{m}(C_k)}$ 

exists, where the limit is taken over all possible partitions ofT of diameter at mostδ.^[3]

Iff is Riemann integrable,S is called theRiemann integral off overT and is denoted

${\displaystyle \int \cdots {\int }_{T}f(x_1,x_2,\dots ,x_n)d{x}_1\cdots d{x}_n}$ .

Frequently this notation is abbreviated as

${\displaystyle {\int }_{T}f(\mathbf{x})d^n\mathbf{x}}$ .

wherex represents then-tuple(x₁, ...,x_n) anddⁿx is then-dimensional volumedifferential.

The Riemann integral of a function defined over an arbitrary boundedn-dimensional set can be defined by extending that function to a function defined over a half-open rectangle whose values are zero outside the domain of the original function. Then the integral of the original function over the original domain is defined to be the integral of the extended function over its rectangular domain, if it exists.

In what follows the Riemann integral inn dimensions will be called themultiple integral.

Multiple integrals have many properties common to those of integrals of functions of one variable (linearity, commutativity, monotonicity, and so on). One important property of multiple integrals is that the value of an integral is independent of the order of integrands under certain conditions. This property is popularly known asFubini's theorem.^[4]

In the case of${\displaystyle T\subseteq {\mathbb{R}}^2}$ , the integral

${\displaystyle l={\iint }_{T}f(x,y)dxdy}$ 

is thedouble integral off onT, and if${\displaystyle T\subseteq {\mathbb{R}}^3}$  the integral

${\displaystyle l={\iiint }_{T}f(x,y,z)dxdydz}$ 

is thetriple integral off onT.

Notice that, by convention, the double integral has two integral signs, and the triple integral has three; this is a notational convention which is convenient when computing a multiple integral as an iterated integral, as shown later in this article.

The resolution of problems with multiple integrals consists, in most cases, of finding a way to reduce the multiple integral to aniterated integral, a series of integrals of one variable, each being directly solvable. For continuous functions, this is justified byFubini's theorem. Sometimes, it is possible to obtain the result of the integration by direct examination without any calculations.

The following are some simple methods of integration:^[1]

When the integrand is aconstant functionc, the integral is equal to the product ofc and the measure of the domain of integration. Ifc = 1 and the domain is a subregion ofR², the integral gives the area of the region, while if the domain is a subregion ofR³, the integral gives the volume of the region.

Example. Letf(x,y) = 2 and

${\displaystyle D=\left\{(x,y)\in {\mathbb{R}}^2:2\le x\le 4;3\le y\le 6\right\}}$ ,

in which case

${\displaystyle {\int }_3^6{\int }_2^{4}2dxdy=2{\int }_3^6{\int }_2^{4}1dxdy=2\cdot \mathrm{area}(D)=2\cdot (2\cdot 3)=12}$ ,

since by definition we have:

${\displaystyle {\int }_3^6{\int }_2^{4}1dxdy=\mathrm{area}(D)}$ .

When the domain of integration is symmetric about the origin with respect to at least one of the variables of integration and the integrand isodd with respect to this variable, the integral is equal to zero, as the integrals over the two halves of the domain have the same absolute value but opposite signs. When the integrand iseven with respect to this variable, the integral is equal to twice the integral over one half of the domain, as the integrals over the two halves of the domain are equal.

Example 1. Consider the functionf(x,y) = 2 sin(x) − 3y³ + 5 integrated over the domain

${\displaystyle T=\left\{(x,y)\in {\mathbb{R}}^2:x^2+y^2\le 1\right\}}$ ,

a disc withradius 1 centered at the origin with the boundary included.

Using the linearity property, the integral can be decomposed into three pieces:

${\displaystyle {\iint }_T\left(2\sin x-3y^3+5\right)dxdy={\iint }_{T}2\sin xdxdy-{\iint }_{T}3y^{3}dxdy+{\iint }_{T}5dxdy}$ .

The function2 sin(x) is an odd function in the variablex and the discT is symmetric with respect to they-axis, so the value of the first integral is 0. Similarly, the function3y³ is an odd function ofy, andT is symmetric with respect to thex-axis, and so the only contribution to the final result is that of the third integral. Therefore the original integral is equal to the area of the disk times 5, or 5π.

Example 2. Consider the functionf(x,y,z) =x exp(y² +z²) and as integration region theball with radius 2 centered at the origin,

${\displaystyle T=\left\{(x,y,z)\in {\mathbb{R}}^3:x^2+y^2+z^2\le 4\right\}}$ .

The "ball" is symmetric about all three axes, but it is sufficient to integrate with respect tox-axis to show that the integral is 0, because the function is an odd function of that variable.

This method is applicable to any domainD for which:

• Theprojection ofD onto either thex-axis or they-axis is bounded by the two values,a andb
• Any line perpendicular to this axis that passes between these two values intersects the domain in an interval whose endpoints are given by the graphs of two functions,α andβ

Such a domain will be here called anormal domain. Elsewhere in the literature, normal domains are sometimes called type I or type II domains, depending on which axis the domain is fibred over. In all cases, the function to be integrated must be Riemann integrable on the domain, which is true (for instance) if the function is continuous.

If the domainD is normal with respect to thex-axis, andf :D →R is acontinuous function; thenα(x) andβ(x) (both of which are defined on the interval[a,b]) are the two functions that determineD. Then, by Fubini's theorem:^[5]

${\displaystyle {\iint }_{D}f(x,y)dxdy={\int }_a^{b}dx{\int }_{\alpha (x)}^{\beta (x)}f(x,y)dy}$ .

IfD is normal with respect to they-axis andf :D →R is a continuous function; thenα(y) andβ(y) (both of which are defined on the interval[a,b]) are the two functions that determineD. Again, by Fubini's theorem:

${\displaystyle {\iint }_{D}f(x,y)dxdy={\int }_a^{b}dy{\int }_{\alpha (y)}^{\beta (y)}f(x,y)dx}$ .

IfT is a domain that is normal with respect to thexy-plane and determined by the functionsα(x,y) andβ(x,y), then

${\displaystyle {\iiint }_{T}f(x,y,z)dxdydz={\iint }_D{\int }_{\alpha (x,y)}^{\beta (x,y)}f(x,y,z)dzdxdy}$ .

This definition is the same for the other five normality cases onR³. It can be generalized in a straightforward way to domains inRⁿ.

The limits of integration are often not easily interchangeable (without normality or with complex formulae to integrate). One makes achange of variables to rewrite the integral in a more "comfortable" region, which can be described in simpler formulae. To do so, the function must be adapted to the new coordinates.

Example 1a. The function isf(x,y) = (x − 1)² +√y; if one adopts the substitutionu =x − 1,v =y thereforex =u + 1,y =v one obtains the new functionf₂(u,v) = (u)² +√v.

• Similarly for the domain because it is delimited by the original variables that were transformed before (x andy in example)
• The differentialsdx anddy transform via the absolute value of thedeterminant of the Jacobian matrix containing the partial derivatives of the transformations regarding the new variable (consider, as an example, the differential transformation in polar coordinates)

There exist three main "kinds" of changes of variable (one inR², two inR³); however, more general substitutions can be made using the same principle.

InR² if the domain has a circular symmetry and the function has some particular characteristics one can apply thetransformation to polar coordinates (see the example in the picture) which means that the generic pointsP(x,y) in Cartesian coordinates switch to their respective points in polar coordinates. That allows one to change the shape of the domain and simplify the operations.

The fundamental relation to make the transformation is the following:

${\displaystyle f(x,y)\to f(\rho \cos \phi ,\rho \sin \phi )}$ .

Example 2a. The function isf(x,y) =x +y and applying the transformation one obtains

${\displaystyle f(x,y)=f(\rho \cos \phi ,\rho \sin \phi )=\rho \cos \phi +\rho \sin \phi =\rho (\cos \phi +\sin \phi )}$ .

Example 2b. The function isf(x,y) =x² +y², in this case one has:

${\displaystyle f(x,y)={\rho }^2\left({\cos }^2\phi +{\sin }^2\phi \right)={\rho }^2}$ 

using thePythagorean trigonometric identity (can be useful to simplify this operation).

The transformation of the domain is made by defining the radius' crown length and the amplitude of the described angle to define theρ,φ intervals starting fromx,y.

Example 2c. The domain isD = {x² +y² ≤ 4}, that is a circumference of radius 2; it's evident that the covered angle is the circle angle, soφ varies from 0 to 2π, while the crown radius varies from 0 to 2 (the crown with the inside radius null is just a circle).

Example 2d. The domain isD = {x² +y² ≤ 9,x² +y² ≥ 4,y ≥ 0}, that is the circular crown in the positivey half-plane (please see the picture in the example);φ describes a plane angle whileρ varies from 2 to 3. Therefore the transformed domain will be the followingrectangle:

${\displaystyle T=\{2\le \rho \le 3,0\le \phi \le \pi \}}$ .

TheJacobian determinant of that transformation is the following:

 ,

which has been obtained by inserting the partial derivatives ofx =ρ cos(φ),y =ρ sin(φ) in the first column respect toρ and in the second respect toφ, so thedx dy differentials in this transformation becomeρ dρ dφ.

Once the function is transformed and the domain evaluated, it is possible to define the formula for the change of variables in polar coordinates:

${\displaystyle {\iint }_{D}f(x,y)dxdy={\iint }_{T}f(\rho \cos \phi ,\rho \sin \phi )\rho d\rho d\phi }$ .

φ is valid in the[0, 2π] interval whileρ, which is a measure of a length, can only have positive values.

Example 2e. The function isf(x,y) =x and the domain is the same as in Example 2d. From the previous analysis ofD we know the intervals ofρ (from 2 to 3) and ofφ (from 0 toπ). Now we change the function:

${\displaystyle f(x,y)=x⟶f(\rho ,\phi )=\rho \cos \phi }$ .

Finally let's apply the integration formula:

${\displaystyle {\iint }_{D}xdxdy={\iint }_T\rho \cos \phi \rho d\rho d\phi }$ .

Once the intervals are known, you have

${\displaystyle {\int }_0^{\pi }{\int }_2^3{\rho }^2\cos \phi d\rho d\phi ={\int }_0^{\pi }\cos \phi d\phi {\left[\frac{{\rho }^3}{3}\right]}_2^3=[\sin \phi {]}_0^{\pi }\left(9-\frac{8}{3}\right)=0}$ .

InR³ the integration on domains with a circular base can be made by thepassage tocylindrical coordinates; the transformation of the function is made by the following relation:

${\displaystyle f(x,y,z)\to f(\rho \cos \phi ,\rho \sin \phi ,z)}$ 

The domain transformation can be graphically attained, because only the shape of the base varies, while the height follows the shape of the starting region.

Example 3a. The region isD = {x² +y² ≤ 9,x² +y² ≥ 4, 0 ≤z ≤ 5} (that is the "tube" whose base is the circular crown of Example 2d and whose height is 5); if the transformation is applied, this region is obtained:

${\displaystyle T=\{2\le \rho \le 3,0\le \phi \le 2\pi ,0\le z\le 5\}}$ 

(that is, the parallelepiped whose base is similar to the rectangle in Example 2d and whose height is 5).

Because thez component is unvaried during the transformation, thedx dy dz differentials vary as in the passage to polar coordinates: therefore, they becomeρ dρ dφ dz.

Finally, it is possible to apply the final formula to cylindrical coordinates:

${\displaystyle {\iiint }_{D}f(x,y,z)dxdydz={\iiint }_{T}f(\rho \cos \phi ,\rho \sin \phi ,z)\rho d\rho d\phi dz}$ .

This method is convenient in case of cylindrical or conical domains or in regions where it is easy to individuate thez interval and even transform the circular base and the function.

Example 3b. The function isf(x,y,z) =x² +y² +z and as integration domain thiscylinder:D = {x² +y² ≤ 9, −5 ≤z ≤ 5}. The transformation ofD in cylindrical coordinates is the following:

${\displaystyle T=\{0\le \rho \le 3,0\le \phi \le 2\pi ,-5\le z\le 5\}}$ .

while the function becomes

${\displaystyle f(\rho \cos \phi ,\rho \sin \phi ,z)={\rho }^2+z}$ .

Finally one can apply the integration formula:

${\displaystyle {\iiint }_D\left(x^2+y^2+z\right)dxdydz={\iiint }_T\left({\rho }^2+z\right)\rho d\rho d\phi dz}$ ;

developing the formula you have

${\displaystyle {\int }_{-5}^{5}dz{\int }_0^{2\pi }d\phi {\int }_0^3\left({\rho }^3+\rho z\right)d\rho =2\pi {\int }_{-5}^5{\left[\frac{{\rho }^4}{4}+\frac{{\rho }^{2}z}{2}\right]}_0^{3}dz=2\pi {\int }_{-5}^5\left(\frac{81}{4}+\frac{9}{2}z\right)dz=\cdots =405\pi }$ .

InR³ some domains have a spherical symmetry, so it's possible to specify the coordinates of every point of the integration region by two angles and one distance. It's possible to use therefore thepassage tospherical coordinates; the function is transformed by this relation:

${\displaystyle f(x,y,z)⟶f(\rho \cos \theta \sin \phi ,\rho \sin \theta \sin \phi ,\rho \cos \phi )}$ .

Points on thez-axis do not have a precise characterization in spherical coordinates, soθ can vary between 0 and 2π.

The better integration domain for this passage is the sphere.

Example 4a. The domain isD =x² +y² +z² ≤ 16 (sphere with radius 4 and center at the origin); applying the transformation you get the region

${\displaystyle T=\{0\le \rho \le 4,0\le \phi \le \pi ,0\le \theta \le 2\pi \}}$ .

The Jacobian determinant of this transformation is the following:

 .

Thedx dy dz differentials therefore are transformed toρ² sin(φ)dρdθdφ.

This yields the final integration formula:

${\displaystyle {\iiint }_{D}f(x,y,z)dxdydz={\iiint }_{T}f(\rho \sin \phi \cos \theta ,\rho \sin \phi \sin \theta ,\rho \cos \phi ){\rho }^2\sin \phi d\rho d\theta d\phi }$ .

It is better to use this method in case of spherical domainsand in case of functions that can be easily simplified by the first fundamental relation of trigonometry extended toR³ (see Example 4b); in other cases it can be better to use cylindrical coordinates (see Example 4c).

${\displaystyle {\iiint }_{T}f(a,b,c){\rho }^2\sin \phi d\rho d\theta d\phi }$ .

The extraρ² andsinφ come from the Jacobian.

In the following examples the roles ofφ andθ have been reversed.

Example 4b.D is the same region as in Example 4a andf(x,y,z) =x² +y² +z² is the function to integrate. Its transformation is very easy:

${\displaystyle f(\rho \sin \phi \cos \theta ,\rho \sin \phi \sin \theta ,\rho \cos \phi )={\rho }^2}$ ,

while we know the intervals of the transformed regionT fromD:

${\displaystyle T=\{0\le \rho \le 4,0\le \phi \le \pi ,0\le \theta \le 2\pi \}}$ .

We therefore apply the integration formula:

${\displaystyle {\iiint }_D\left(x^2+y^2+z^2\right)dxdydz={\iiint }_T{\rho }^2{\rho }^2\sin \theta d\rho d\theta d\phi }$ ,

and, developing, we get

${\displaystyle {\iiint }_T{\rho }^4\sin \theta d\rho d\theta d\phi ={\int }_0^{\pi }\sin \phi d\phi {\int }_0^4{\rho }^{4}d\rho {\int }_0^{2\pi }d\theta =2\pi {\int }_0^{\pi }\sin \phi {\left[\frac{{\rho }^5}{5}\right]}_0^{4}d\phi =2\pi {\left[\frac{{\rho }^5}{5}\right]}_0^4[-\cos \phi {]}_0^{\pi }=\frac{4096\pi }{5}}$ .

Example 4c. The domainD is the ball with center at the origin and radius3a,

${\displaystyle D=\left\{x^2+y^2+z^2\le 9a^2\right\}}$ ,

andf(x,y,z) =x² +y² is the function to integrate.

Looking at the domain, it seems convenient to adopt the passage to spherical coordinates, in fact, the intervals of the variables that delimit the newT region are:

${\displaystyle T=\{0\le \rho \le 3a,0\le \phi \le 2\pi ,0\le \theta \le \pi \}}$ .

However, applying the transformation, we get

${\displaystyle f(x,y,z)=x^2+y^2⟶{\rho }^2{\sin }^2\theta {\cos }^2\phi +{\rho }^2{\sin }^2\theta {\sin }^2\phi ={\rho }^2{\sin }^2\theta }$ .

Applying the formula for integration we obtain:

${\displaystyle {\iiint }_T{\rho }^2{\sin }^2\theta {\rho }^2\sin \theta d\rho d\theta d\phi ={\iiint }_T{\rho }^4{\sin }^3\theta d\rho d\theta d\phi }$ ,

which can be solved by turning it into an iterated integral.${\displaystyle {\iiint }_T{\rho }^4{\sin }^3\theta d\rho d\theta d\phi =\underset{I}{\underbrace{{\int }_0^{3a}{\rho }^{4}d\rho }}\underset{II}{\underbrace{{\int }_0^{\pi }{\sin }^3\theta d\theta }}\underset{III}{\underbrace{{\int }_0^{2\pi }d\phi }}}$ .

${\displaystyle I={{\int }_0^{3a}{\rho }^{4}d\rho =\frac{{\rho }^5}{5}|}_0^{3a}=\frac{243}{5}{a}^5}$ ,

${\displaystyle II={\int }_0^{\pi }{\sin }^3\theta d\theta =-{\int }_0^{\pi }{\sin }^2\theta d(\cos \theta )={\int }_0^{\pi }({\cos }^2\theta -1)d(\cos \theta )={\frac{{\cos }^3\theta }{3}|}_0^{\pi }-{\cos \theta |}_0^{\pi }=\frac{4}{3}}$ ,

${\displaystyle III={\int }_0^{2\pi }d\phi =2\pi }$ .

Collecting all parts,

${\displaystyle {\iiint }_T{\rho }^4{\sin }^3\theta d\rho d\theta d\phi =I\cdot II\cdot III=\frac{243}{5}{a}^5\cdot \frac{4}{3}\cdot 2\pi =\frac{648}{5}\pi a^5}$ .

Alternatively, this problem can be solved by using the passage to cylindrical coordinates. The newT intervals are

${\displaystyle T=\left\{0\le \rho \le 3a,0\le \phi \le 2\pi ,-\sqrt{9a^2-{\rho }^2}\le z\le \sqrt{9a^2-{\rho }^2}\right\}}$ ;

thez interval has been obtained by dividing the ball into twohemispheres simply by solving theinequality from the formula ofD (and then directly transformingx² +y² intoρ²). The new function is simplyρ². Applying the integration formula

${\displaystyle {\iiint }_T{\rho }^2\rho d\rho d\phi dz}$ .

Then we get:

 

Thanks to the passage to cylindrical coordinates it was possible to reduce the triple integral to an easier one-variable integral.

See also the differential volume entry innabla in cylindrical and spherical coordinates.

Let us assume that we wish to integrate a multivariable functionf over a regionA:

${\displaystyle A=\left\{(x,y)\in {\mathbf{R}}^2:11\le x\le 14;7\le y\le 10\right\}\text{ and }f(x,y)=x^2+4y}$ .

From this we formulate the iterated integral

${\displaystyle {\int }_7^{10}{\int }_{11}^{14}(x^2+4y)dxdy}$ .

The inner integral is performed first, integrating with respect tox and takingy as a constant, as it is not thevariable of integration. The result of this integral, which is a function depending only ony, is then integrated with respect toy.

 

We then integrate the result with respect toy.

 

In cases where the double integral of the absolute value of the function is finite, the order of integration is interchangeable, that is, integrating with respect tox first and integrating with respect toy first produce the same result. That isFubini's theorem. For example, doing the previous calculation with order reversed gives the same result:

 

Consider the region (please see the graphic in the example):

${\displaystyle D=\{(x,y)\in {\mathbf{R}}^2:x\ge 0,y\le 1,y\ge x^2\}}$  .

Calculate

${\displaystyle {\iint }_D(x+y)dxdy}$ .

This domain is normal with respect to both thex- andy-axes. To apply the formulae it is required to find the functions that determineD and the intervals over which these functions are defined. In this case the two functions are:

${\displaystyle \alpha (x)=x^2\text{ and }\beta (x)=1}$ 

while the interval is given by the intersections of the functions withx = 0, so the interval is [a, b] = [0, 1] (normality has been chosen with respect to thex-axis for a better visual understanding).

It is now possible to apply the formula:

${\displaystyle {\iint }_D(x+y)dxdy={\int }_0^{1}dx{\int }_{x^2}^1(x+y)dy={\int }_0^{1}dx{\left[xy+\frac{y^2}{2}\right]}_{x^2}^1}$ 

(at first the second integral is calculated consideringx as a constant). The remaining operations consist of applying the basic techniques of integration:

${\displaystyle {\int }_0^1{\left[xy+\frac{y^2}{2}\right]}_{x^2}^{1}dx={\int }_0^1\left(x+\frac{1}{2}-x^3-\frac{x^4}{2}\right)dx=\cdots =\frac{13}{20}}$ .

If we choose normality with respect to they-axis we could calculate

${\displaystyle {\int }_0^{1}dy{\int }_0^{\sqrt{y}}(x+y)dx}$ .

and obtain the same value.

Using the methods previously described, it is possible to calculate the volumes of some common solids.

• Cylinder: The volume of a cylinder with heighth and circular base of radiusR can be calculated by integrating the constant functionh over the circular base, using polar coordinates.

${\displaystyle \mathrm{V}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{e}={\int }_0^{2\pi }d\phi {\int }_0^{R}h\rho d\rho =2\pi h{\left[\frac{{\rho }^2}{2}\right]}_0^R=\pi R^{2}h}$ 

This is in agreement with the formula for the volume of aprism

${\displaystyle \mathrm{V}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{e}=\text{base area}\times \text{height}}$ .

• Sphere: The volume of a sphere with radiusR can be calculated by integrating the constant function 1 over the sphere, using spherical coordinates.

 

• Tetrahedron (triangularpyramid or 3-simplex): The volume of a tetrahedron with its apex at the origin and edges of lengthℓ along thex-,y- andz-axes can be calculated by integrating the constant function 1 over the tetrahedron.

 

This is in agreement with the formula for the volume of apyramid.

${\displaystyle \mathrm{V}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{e}=\frac{1}{3}\times \text{base area}\times \text{height}=\frac{1}{3}\times \frac{{\ell }^2}{2}\times \ell =\frac{{\ell }^3}{6}}$ .

In case of unbounded domains or functions not bounded near the boundary of the domain, we have to introduce thedoubleimproper integral or thetriple improper integral.

Fubini's theorem states that if^[4]

 ,

that is, if the integral is absolutely convergent, then the multiple integral will give the same result as either of the two iterated integrals:

${\displaystyle {\iint }_{A\times B}f(x,y)d(x,y)={\int }_A\left({\int }_{B}f(x,y)dy\right)dx={\int }_B\left({\int }_{A}f(x,y)dx\right)dy}$ .

In particular this will occur if|f(x,y)| is abounded function andA andB arebounded sets.

If the integral is not absolutely convergent, care is needed not to confuse the concepts ofmultiple integral anditerated integral, especially since the same notation is often used for either concept. The notation

${\displaystyle {\int }_0^1{\int }_0^{1}f(x,y)dydx}$ 

means, in some cases, an iterated integral rather than a true double integral. In an iterated integral, the outer integral

${\displaystyle {\int }_0^1\cdots dx}$ 

is the integral with respect tox of the following function ofx:

${\displaystyle g(x)={\int }_0^{1}f(x,y)dy}$ .

A double integral, on the other hand, is defined with respect to area in thexy-plane. If the double integral exists, then it is equal to each of the two iterated integrals (either "dy dx" or "dx dy") and one often computes it by computing either of the iterated integrals. But sometimes the two iterated integrals exist when the double integral does not, and in some such cases the two iterated integrals are different numbers, i.e., one has

${\displaystyle {\int }_0^1{\int }_0^{1}f(x,y)dydx\ne {\int }_0^1{\int }_0^{1}f(x,y)dxdy}$ .

This is an instance of rearrangement of aconditionally convergent integral.

On the other hand, some conditions ensure that the two iterated integrals are equal even though the double integral need not exist. By theFichtenholz–Lichtenstein theorem, iff is bounded on[0, 1] × [0, 1] and both iterated integrals exist, then they are equal. Moreover, existence of the inner integrals ensures existence of the outer integrals.^[6][7][8] The double integral need not exist in this case even asLebesgue integral, according toSierpiński.^[9]

The notation

${\displaystyle {\int }_{[0,1]\times [0,1]}f(x,y)dxdy}$ 

may be used if one wishes to be emphatic about intending a double integral rather than an iterated integral.

Triple integral was demonstrated by Fubini's theorem.^[10] Drichlet theorem and Liouville 's extension theorem on Triple integral.

Quite generally, just as in one variable, one can use the multiple integral to find the average of a function over a given set. Given a setD ⊆Rⁿ and an integrable functionf overD, the average value off over its domain is given by

${\displaystyle \overline{f}=\frac{1}{m(D)}{\int }_{D}f(x)dx}$ ,

wherem(D) is themeasure ofD.

Additionally, multiple integrals are used in many applications inphysics. The examples below also show some variations in the notation.

Inmechanics, themoment of inertia is calculated as the volume integral (triple integral) of thedensity weighed with the square of the distance from the axis:

${\displaystyle I_z={\iiint }_V\rho r^{2}dV}$ .

Thegravitational potential associated with amass distribution given by a massmeasuredm on three-dimensionalEuclidean spaceR³ is^[11]

${\displaystyle V(\mathbf{x})=-{\iiint }_{{\mathbf{R}}^3}\frac{G}{|\mathbf{x}-\mathbf{y}|}dm(\mathbf{y})}$ .

If there is a continuous functionρ(x) representing the density of the distribution atx, so thatdm(x) =ρ(x)d³x, whered³x is the Euclideanvolume element, then the gravitational potential is

${\displaystyle V(\mathbf{x})=-{\iiint }_{{\mathbf{R}}^3}\frac{G}{|\mathbf{x}-\mathbf{y}|}\rho (\mathbf{y})d^3\mathbf{y}}$ .

Inelectromagnetism,Maxwell's equations can be written using multiple integrals to calculate the total magnetic and electric fields.^[12] In the following example, theelectric field produced by a distribution ofcharges given by the volumecharge densityρ(r→ ) is obtained by atriple integral of a vector function:

${\displaystyle \overrightarrow{E}=\frac{1}{4\pi {\epsilon }_0}\iiint \frac{\overrightarrow{r}-{\overrightarrow{r}}^{\prime }}{{\Vert \overrightarrow{r}-{\overrightarrow{r}}^{\prime }\Vert }^3}\rho ({\overrightarrow{r}}^{\prime })d^3{r}^{\prime }}$ .

This can also be written as an integral with respect to asigned measure representing the charge distribution.

• Mainanalysis theorems that relate multiple integrals:
  • Divergence theorem
  • Stokes' theorem
  • Green's theorem

• Adams, Robert A. (2003).Calculus: A Complete Course (5th ed.). Addison-Wesley.ISBN 0-201-79131-5.
• Jain, R. K.; Iyengar, S. R. K. (2009).Advanced Engineering Mathematics (3rd ed.). Narosa Publishing House.ISBN 978-81-7319-730-7.
• Herman, Edwin “Jed” & Strang, Gilbert (2016):Calculus : Volume 3 : OpenStax, Rice University, Houston, Texas, USA.ISBN 978-1-50669-805-2. (PDF)

• Weisstein, Eric W."Multiple Integral".MathWorld.
• L.D. Kudryavtsev (2001) [1994],"Multiple integral",Encyclopedia of Mathematics,EMS Press
• Mathematical Assistant on Web online evaluation of double integrals inCartesian coordinates andpolar coordinates (includes intermediate steps in the solution, powered byMaxima (software))
• Online Double Integral Calculator by WolframAlpha
• Online Triple Integral Calculator by WolframAlpha