Double Integral Calculator

A double integral ∫∫_R f(x,y) dA computes the volume under the surface z = f(x,y) above a region R in the xy-plane. For a rectangular region [a,b] × [c,d], it can be evaluated as an iterated integral:

∫_a^b ∫_c^d f(x,y) dy dx

Inner integral first (treating x as constant), then outer integral.

This calculator uses Simpson’s rule in two dimensions. For each of n+1 evenly spaced x values, it integrates f(x,y) over y using 1D Simpson’s rule, giving a function g(x). Then it integrates g(x) over the x range, again with Simpson’s rule. The combined error is O(h⁴) in the step size h, which at n = 40 puts a smooth integrand right at the edge of double precision.

All three integrand forms here also have a closed-form antiderivative, so the result panel prints the exact value beside the numerical one. If the two disagree by more than rounding, that is worth knowing, and it is the reason both are shown rather than just the answer.

Three integrand forms are available:

Sum: f(x,y) = ax^p + by^q. Easily verifiable by hand for small p, q. Product: f(x,y) = ax^p * y^q. Factors into ∫x^p dx · ∫y^q dy when integrated over a rectangle. Constant: f(x,y) = a. Integral equals a × (area of rectangle).

Worked example: ∫_0^2 ∫_0^1 (x² + y²) dy dx. Inner: ∫_0^1 (x² + y²) dy = x²·y + y³/3 from 0 to 1 = x² + 1/3. Outer: ∫_0^2 (x² + 1/3) dx = x³/3 + x/3 from 0 to 2 = 8/3 + 2/3 = 10/3 ≈ 3.333.

Use cases: volume under a surface, mass with variable density, average value of a function over a region, center of mass coordinates, moments of inertia. In probability, double integrals compute probabilities for joint continuous distributions.

Non-rectangular regions need variable bounds (Type I or Type II regions), or a change to polar or elliptic coordinates with a Jacobian. Neither is in scope for a fixed-rectangle calculator, so if your region is a disc or a triangle this is the wrong tool.