Multiple integration

Multivariable Calculus, Math 211, Spring 2017, Carleton College, Prof. Joshua R. Davis

Infrastructure for Riemann sums

Activate this code cell, which defines functions for computing and drawing the Riemann sum that approximates a Riemann integral. You don’t have to understand this code.

Clear[x, y, boxes, integrate];
boxes[f_, a_, b_, c_, d_, n_] :=
	Module[{dx, dy, max, min, polygons, lines, i, j, xi, yj, xi1, yj1, fij},
		dx = (b - a) / n;
		dy = (d - c) / n;
		min = f /. {x -> b, y -> d};
		max = min;
		polygons = {};
		lines = {};
		For[i = 1, i ≤ n, i += 1,
			xi = a + i dx;
			xi1 = a + (i - 1) dx;
			For[j = 1, j ≤ n, j += 1,
				yj = c + j dy;
				yj1 = c + (j - 1) dy;
				fij = f /. {x -> xi, y -> yj};
				min = Min[min, fij];
				max = Max[max, fij];
				polygons = Append[polygons,
					{{xi1, yj1, fij}, {xi, yj1, fij}, {xi, yj, fij}, {xi1, yj, fij}}];
				lines = Join[lines, {
					{{xi, yj, 0}, {xi, yj, fij}},
					{{xi1, yj, 0}, {xi1, yj, fij}},
					{{xi1, yj1, 0}, {xi1, yj1, fij}},
					{{xi, yj1, 0}, {xi, yj1, fij}}}]]];
		Show[
			Graphics3D[ Join[Map[Line, lines], Map[Polygon, polygons]], AspectRatio-> 1],
			Plot3D[f, {x, a, b}, {y, c, d}, PlotRange->{min, max}]]];
integrate[f_, a_, b_, c_, d_, n_] :=
	Module[{dx, dy, sum, i, j, xi, yj, xi1, yj1, fij},
		dx = (b - a) / n;
		dy = (d - c) / n;
		sum = 0;
		For[i = 1, i ≤ n, i += 1,
			xi = a + i dx;
			For[j = 1, j ≤ n, j += 1,
				yj = c + j dy;
				fij = f /. {x -> xi, y -> yj};
				sum += fij dx dy]];
		sum];

Plane

Try changing n, to see how the approximated integral changes. Gradually it approaches the true value, which Mathematica can compute for us (using a more advanced version of the same idea).

Clear[x, y];
f = 2 x + y + 1;
a = -1;
b = 1;
c = -1;
d = 1;
n = 16;
N[integrate[f, a, b, c, d, n]]
NIntegrate[f, {x, a, b}, {y, c, d}]
boxes[f, a, b, c, d, n]

Monkey saddle

Try changing n, to see how the approximated integral changes. Gradually it approaches the true value of 0. How can you tell, just by looking at the function, that the integral equals 0?

Clear[x, y];
f = 3 x y^2 - x^3;
a = -1;
b = 1;
c = -1;
d = 1;
n = 8;
N[integrate[f, a, b, c, d, n]]
NIntegrate[f, {x, a, b}, {y, c, d}]
boxes[f, a, b, c, d, n]

Integration in practice

The following two commands both compute a double integral of f(x, y) over a certain region, using iterated integrals. The first command computes the iterated integrals in dx dy order, and the second computes them in dy dx order. In view of Fubini’s theorem, does the answer make sense?

Integrate[x Sin[y], {x, 0, 1}, {y, -1, 3}]
Integrate[x Sin[y], {y, -1, 3}, {x, 0, 1}]

Here’s another example. In view of Fubini’s theorem, does the answer make sense?

Integrate[(x - y) / (x^3 + y^3), {x, 0, 1}, {y, 0, 1}]
Integrate[(x - y) / (x^3 + y^3), {y, 0, 1}, {x, 0, 1}]

Mimicking the “integrate” examples above, compute a Riemann sum approximation to this double integral, to get another idea of what the answer is.