What Is Numerical Evaluation of Double Integrals?
A double integral measures accumulation over a two-dimensional region. Instead of adding values along a line, as a single integral does, a double integral adds values over an area in the \(xy\)-plane.
For a function \(f(x,y)\) over a region \(R\), the double integral is written as:
$$
\iint_R f(x,y)\,dA
$$
When \(f(x,y)\) is nonnegative and represents height, this integral can be interpreted as the volume under the surface \(z=f(x,y)\) and above the region \(R\). When the function has both positive and negative values, the result is a signed accumulation: positive values add to the total, while negative values subtract from it.
Many double integrals can be evaluated exactly by hand, but not all of them are convenient or even possible to express using elementary antiderivatives. Numerical evaluation solves this practical problem by estimating the integral from a finite set of sample points.
This calculator begins with midpoint sampling. It then doubles the grid and uses Richardson extrapolation to compare successive estimates. A result is shown only when the adaptive estimates agree within the calculator's numerical tolerance.
Why Double Integrals Matter
Double integrals are useful because many quantities are spread across a region rather than located at one point. In calculus, physics, engineering, statistics, and applied modeling, they can represent ideas such as:
-
volume under a surface;
-
total mass from a density spread across a plate;
-
total charge, heat, rainfall, or population density over a region;
-
average value of a function across an area;
-
signed accumulation when values above and below zero cancel.
Numerical double integration is especially useful when the region is easy to describe but the antiderivative is difficult, when the function is complicated, or when a student wants to check whether an integral setup is reasonable.
Key Terms to Know
-
Integrand: The function being integrated, such as \(f(x,y)\).
-
Region of integration: The area in the \(xy\)-plane over which values are accumulated.
-
Iterated integral: A double integral written as two single-variable integrations, one inside the other.
-
Riemann sum: A finite sum that approximates an integral by multiplying sample values by small areas.
-
Midpoint rule: A numerical method that samples each small cell at its center.
-
Polar coordinates: A coordinate system using radius \(r\) and angle \(\theta\) instead of \(x\) and \(y\).
-
Jacobian factor: The area-scaling factor needed when changing coordinates. In polar coordinates, this factor is \(r\).
-
Average value: The integral divided by the area of the region.
How Double Integrals and Midpoint Sampling Work
For a rectangular region
$$
R=[a,b]\times[c,d],
$$
the exact double integral can be written as:
$$
\int_a^b \int_c^d f(x,y)\,dy\,dx
$$
A midpoint approximation divides the \(x\)-interval and \(y\)-interval into \(n\) parts:
$$
\Delta x=\frac{b-a}{n}, \qquad \Delta y=\frac{d-c}{n}
$$
The midpoint sample coordinates are:
$$
x_i^*=a+\left(i-\frac{1}{2}\right)\Delta x, \qquad y_j^*=c+\left(j-\frac{1}{2}\right)\Delta y
$$
The rectangular midpoint approximation is:
$$
\iint_R f(x,y)\,dA \approx \sum_{i=1}^{n}\sum_{j=1}^{n} f(x_i^*,y_j^*)\,\Delta x\,\Delta y
$$
This means each sample value is treated as the representative height for one small rectangle.
Adaptive refinement and error estimation
Let \(M_n\) be the midpoint estimate from an \(n\)-by-\(n\) grid. Because the midpoint rule is second order for sufficiently smooth inputs, the calculator combines two grid levels using Richardson extrapolation:
$$
R_{2n}=M_{2n}+\frac{M_{2n}-M_n}{3}
$$
It then compares consecutive extrapolated values. The displayed internal error estimate is based on \(\lvert R_{4n}-R_{2n}\rvert/15\). This check catches many unstable, singular, and undersampled inputs, but it is not a rigorous bound for every possible function.
Custom regions with bounds depending on \(x\)
A common nonrectangular region is described by \(x\) running from \(a\) to \(b\), while \(y\) runs from a lower curve \(g(x)\) to an upper curve \(h(x)\):
$$
\int_a^b \int_{g(x)}^{h(x)} f(x,y)\,dy\,dx
$$
For this type of region, each vertical slice can have a different height. The calculator samples an \(x\) midpoint, computes the lower and upper \(y\)-bounds there, then applies midpoint sampling within that slice.
At an \(x\) midpoint \(x_i^*\), the vertical cell height is:
$$
\Delta y_i=\frac{h(x_i^*)-g(x_i^*)}{n}
$$
A corresponding midpoint sample in that slice is:
$$
y_{ij}^*=g(x_i^*)+\left(j-\frac{1}{2}\right)\Delta y_i
$$
The custom-region approximation follows the same idea:
$$
\sum_{i=1}^{n}\sum_{j=1}^{n} f(x_i^*,y_{ij}^*)\,\Delta x\,\Delta y_i
$$
Polar regions and the extra factor \(r\)
Polar coordinates are often simpler for disks, circles, annular regions, sectors, and regions with circular symmetry. The coordinate conversion is:
$$
x=r\cos(\theta), \qquad y=r\sin(\theta)
$$
In polar form, area is not just \(dr\,d\theta\). The area element is:
$$
dA=r\,dr\,d\theta
$$
So a polar double integral is written as:
$$
\int_{\theta=c}^{d}\int_{r=a}^{b} f(r\cos\theta,r\sin\theta)\,r\,dr\,d\theta
$$
The midpoint approximation becomes:
$$
\sum_{i=1}^{n}\sum_{j=1}^{n} f(r_i^*\cos\theta_j^*,r_i^*\sin\theta_j^*)\,r_i^*\,\Delta r\,\Delta\theta
$$
The factor \(r_i^*\) is important. It accounts for the fact that polar cells farther from the origin cover more area than polar cells near the origin.
Area and average value
The area of a region can be found by integrating \(1\) over that region:
$$
A(R)=\iint_R 1\,dA
$$
The average value of \(f\) over a region of positive area is:
$$
f_{\text{avg}}=\frac{1}{A(R)}\iint_R f(x,y)\,dA
$$
The average value has the same units as the function. The double integral has function-value times area units.
Examples of Double Integrals in Practice
Example 1: Rectangular midpoint sampling
Suppose:
$$
f(x,y)=x+y
$$
on the rectangle
$$
0\le x\le 2, \qquad 0\le y\le 1
$$
Using \(n=2\) samples per dimension gives:
$$
\Delta x=1, \qquad \Delta y=0.5, \qquad \Delta A=0.5
$$
The four midpoint values are:
$$
f(0.5,0.25)=0.75
$$
$$
f(0.5,0.75)=1.25
$$
$$
f(1.5,0.25)=1.75
$$
$$
f(1.5,0.75)=2.25
$$
The midpoint estimate is:
$$
(0.75+1.25+1.75+2.25)(0.5)=3
$$
For this linear function, the midpoint estimate matches the exact integral:
$$
\int_0^2\int_0^1 (x+y)\,dy\,dx=3
$$
Example 2: Polar integral over a disk
For the unit disk, use:
$$
0\le r\le 1, \qquad 0\le \theta\le 2\pi
$$
Let:
$$
f(x,y)=x^2+y^2
$$
Since \(x^2+y^2=r^2\), the polar integral is:
$$
\int_0^{2\pi}\int_0^1 r^2\cdot r\,dr\,d\theta
$$
The second \(r\) is the polar area factor, not part of the original function. Evaluating gives:
$$
\int_0^{2\pi}\left[\frac{r^4}{4}\right]_0^1 d\theta
=\int_0^{2\pi}\frac{1}{4}\,d\theta
=\frac{\pi}{2}
$$
In polar mode, the function should still be entered in terms of \(x\) and \(y\). The coordinate conversion and the polar area factor are handled as part of the polar setup.
Example 3: Custom region between two curves
Suppose the region lies between:
$$
y=x^2 \quad \text{and} \quad y=x
$$
from \(x=0\) to \(x=1\). If \(f(x,y)=1\), then the double integral gives the area of the region:
$$
\int_0^1\int_{x^2}^{x}1\,dy\,dx
$$
Evaluate the inner integral first:
$$
\int_{x^2}^{x}1\,dy=x-x^2
$$
Then integrate with respect to \(x\):
$$
\int_0^1(x-x^2)\,dx=\left[\frac{x^2}{2}-\frac{x^3}{3}\right]_0^1=\frac{1}{6}
$$
This illustrates how custom \(y\)-bounds describe a region whose vertical height changes with \(x\).
How to Interpret the Result
The main numerical result is an approximation of the signed double integral over the selected region. A positive result means positive sampled values dominate the accumulation. A negative result means negative sampled values dominate. A result near zero can mean the function is small, the region is small, or positive and negative contributions mostly cancel.
The region area is the area of the integration domain in the \(xy\)-plane. It is not the surface area of \(z=f(x,y)\).
The average value is:
$$
\frac{\text{double integral}}{\text{region area}}
$$
It represents the mean value of \(f(x,y)\) over the region. If the function represents height, the average value can be thought of as the constant height that would give the same total volume over the same base region.
The sampled value range shows the smallest and largest function values found at the sample points. This is useful for a quick check, but it is not a proof of the true global minimum or maximum over the entire region.
For polar regions, the result already includes the polar factor \(r\). Do not multiply the function by \(r\) yourself unless the mathematical integrand truly contains that factor in addition to the polar area element.
Common Mistakes and Misconceptions
Entering \(r\) or \(\theta\) in the function field: The function is entered as \(f(x,y)\). In polar mode, the calculator evaluates that function at \(x=r\cos(\theta)\) and \(y=r\sin(\theta)\).
Multiplying by \(r\) manually in polar mode: The polar area element already includes \(r\). Manually adding another \(r\) changes the integral.
Confusing region area with surface area: The region area is the flat area in the input domain. Surface area of \(z=f(x,y)\) is a different calculation.
Expecting an exact symbolic answer: Midpoint sampling gives a numerical approximation. It does not find antiderivatives or simplify symbolic expressions.
Reversing lower and upper bounds: Lower bounds must be less than upper bounds. In custom mode, the upper function \(h(x)\) should be above the lower function \(g(x)\) over a positive-area part of the interval.
Using unsupported variables or notation: The function uses \(x\) and \(y\). Decimal values should use a decimal point, not a decimal comma.
Assuming the sampled value range is exact: The displayed minimum and maximum are based on sampled points only. A function may reach a higher or lower value between sample points.
When to Use Numerical Double Integration
Use numerical double integration when you need to estimate accumulation over a two-dimensional region and an exact symbolic calculation is inconvenient, unnecessary, or hard to obtain.
It is especially helpful for:
-
checking a double-integral setup for calculus homework or instruction;
-
estimating volume under a surface;
-
finding an approximate average value over a region;
-
working with rectangular regions, circular or sector-shaped polar regions, or regions between two \(y\)-bounds;
-
visualizing how the selected region and surface relate to the integral;
-
reviewing an automatic convergence check and estimated numerical error.
For formal coursework, engineering, physics, or other high-stakes applications, use numerical estimates as part of a broader checking process. Compare with exact methods, error estimates, or appropriate professional tools when accuracy requirements matter.
Limitations and Things to Keep in Mind
This calculator uses adaptive composite midpoint refinement with Richardson extrapolation. It does not compute exact symbolic integrals, and its reported numerical error is an internal estimate rather than a guaranteed mathematical error bound.
The entered sample density is the starting grid. The calculator doubles it automatically, up to \(480\) samples per dimension, and suppresses the main result if consecutive extrapolated estimates do not converge. Functions with sharp spikes, discontinuities, singularities, rapid oscillation, or narrow important features can still evade finite sampling and require a specialized method.
The function must evaluate to finite numbers at the sampled points. Expressions that produce non-finite values such as undefined results or infinity are not treated as improper integrals.
Rectangular and polar bounds are constant bounds. Custom mode supports regions of the form:
$$
a\le x\le b, \qquad g(x)\le y\le h(x)
$$
It does not represent every possible irregular region, multiple disconnected \(y\)-intervals for the same \(x\), or regions better described with \(x\) as a function of \(y\).
In custom mode, the upper curve must remain at or above the lower curve. The calculator checks this at the interval endpoints and at its finite set of sampled \(x\)-positions; any representably distinct finite upper value below the lower value at a checked position is rejected instead of silently clipped or treated as an oriented integral. A narrow crossing between checked positions can evade any finite sampling grid, so analyze unusual or closely spaced bounds independently.
Polar radius lower bounds must be nonnegative, and \(\theta\) values are interpreted in radians.
Starting sample density must be a whole number from \(12\) through \(120\). If the field is left blank, the calculator starts at \(48\). Invalid, fractional, or out-of-range values are rejected rather than silently changed.
Small nonzero results are preserved and may be shown in scientific notation. The calculator displays \(0\) only when the computed floating-point result is zero. Inputs whose area, samples, result, or midpoint coordinates exceed the supported floating-point resolution are rejected explicitly; shift or rescale large-offset bounds before retrying.
Graph previews are visual aids. They help you inspect the region and surface, but they are not a substitute for checking the bounds, function, and numerical accuracy.
How to Use This Calculator
-
Enter the function in terms of \(x\) and \(y\).
-
Choose the integration region type: rectangular, polar, or custom.
-
Enter the limits shown for the selected region type.
-
For custom mode, enter \(g(x)\) as the lower \(y\)-bound and \(h(x)\) as the upper \(y\)-bound.
-
Set the starting samples per dimension, or leave the field blank to start at \(48\). The calculator refines the grid automatically.
-
Review the numerical result, region area, average value, sampled value range, estimated numerical error, and final grid density.
-
Use the formula and interpretation area to check the integral setup.
-
Switch between the 2D region and 3D surface views to visually inspect the domain and function behavior.
The function field supports \(x\), \(y\), common constants such as \(\pi\) and \(e\), and common functions such as trigonometric, logarithmic, exponential, absolute-value, square-root, minimum, and maximum functions.
Frequently Asked Questions
Does this calculator give an exact double integral?
No. It gives a midpoint-rule numerical approximation. For exact symbolic work, you still need algebraic integration methods, a computer algebra system, or a formal solution process.
Why do I enter \(f(x,y)\) even in polar mode?
The function field is based on Cartesian variables. In polar mode, the calculator substitutes \(x=r\cos(\theta)\) and \(y=r\sin(\theta)\) when it samples the function.
Should I include the polar factor \(r\) in my function?
No, not when \(r\) is only the polar area factor. Polar mode already multiplies by \(r\) as part of \(dA=r\,dr\,d\theta\). Only include an extra factor if the actual function you want to integrate contains that factor for a separate mathematical reason.
What does the average value mean?
The average value is the double integral divided by the region area. It is the constant function value that would produce the same total accumulation over the same region.
Why is my result negative or close to zero?
A double integral is signed. Negative function values subtract from the total, so positive and negative regions can cancel. A near-zero result does not always mean the function is zero everywhere.
Does the value range show the true minimum and maximum?
Not necessarily. It shows the minimum and maximum values found among the sampled points. The true extrema could occur between sample points or on a boundary that the midpoint grid does not sample directly.
Sources and References
Books and Open Textbooks
-
Gilbert Strang and Edwin “Jed” Herman. Calculus Volume 3. OpenStax, 2016. Chapter 5, especially 5.1 Double Integrals over Rectangular Regions, 5.2 Double Integrals over General Regions, and 5.3 Double Integrals in Polar Coordinates.
-
Matthew Boelkins, David Austin, and Steven Schlicker. Active Calculus Multivariable. Grand Valley State University Libraries, 2017; LibreTexts edition. Chapter 11, especially 11.3 Double Integrals over General Regions and 11.5 Double Integrals in Polar Coordinates.
-
Michael Corral. Vector Calculus. Schoolcraft College, 2008; latest online version listed August 15, 2022. Chapter 3, Multiple Integrals and Numerical Approximation of Multiple Integrals. Book homepage and LibreTexts edition.
Online Course Notes
-
MIT OpenCourseWare. “Part A: Double Integrals.” 18.02SC Multivariable Calculus, accessed June 28, 2026. Course page.
-
MIT OpenCourseWare. “Changing Variables in Multiple Integrals.” 18.02 Multivariable Calculus, Fall 2007 notes, accessed June 28, 2026. PDF notes.
