Seite - 79 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python

3.7 DoubleandTriple Integrals 79 The expression looks somewhat different from (3.21), but that is because of the notation: sincewe integrate in they direction andwill have toworkwith bothx and y as coordinates, wemust use ny for n, hy for h, and the counter i is more naturallycalledjwhenintegratinginy. Integrals in thexdirectionwillusehx and nx forhandn, and i as counter. Thedouble integral is Rb a g.x/dx, which canbe approximatedby themidpoint method: bZ a g.x/dx hx nx 1X iD0 g.xi/; xi DaC 1 2 hx C ihx : Puttingtheformulastogether,wearriveatthecompositemidpointmethodforadou- ble integral: bZ a dZ c f.x;y/dydx hx nx 1X iD0 hy ny 1X jD0 f.xi;yj/ Dhxhy nx 1X iD0 ny 1X jD0 f aC hx 2 C ihx;cC hy 2 Cjhy : (3.25) Direct derivation The formula (3.25) can also be derived directly in the two- dimensional case by applying the idea of the midpoint method. We divide the rectangle Œa;b Œc;d intonx ny equal-sized cells. The idea of themidpoint method is toapproximatef byaconstantovereachcell, andevaluate theconstant at themidpoint. Cell .i;j/occupies thearea ŒaC ihx;aC.iC1/hx ŒcCjhy;cC.j C1/hy ; and themidpoint is .xi;yj/with xi DaC ihx C 1 2 hx; yj DcCjhy C 1 2 hy : The integral over the cell is thereforehxhyf.xi;yj/, and the total double integral is thesumoverall cells,which isnothingbut formula (3.25). Programming a double sum The formula (3.25) involves a double sum, which is normally implemented as a double for loop. A Python function implementing (3.25)may look like def midpoint_double1(f, a, b, c, d, nx, ny): hx = (b - a)/float(nx) hy = (d - c)/float(ny) I = 0 for i in range(nx): for j in range(ny): xi = a + hx/2 + i*hx yj = c + hy/2 + j*hy I += hx*hy*f(xi, yj) return I