Page - 157 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python 3.6, Volume Second Edition

6.7 DoubleandTriple Integrals 157 Demonstrating Correct Convergence Rates Computing convergence rates re- quires somewhat more tedious programming than for the previous tests, but it can beapplied tomoregeneral integrands.Thealgorithmtypicallygoes like • for i=0,1,2,.. .,q – ni =2i+1 – Compute integralwithni intervals – Compute theerrorEi – Estimate ri from(6.28) if i >0 Thecorrespondingcodemaylook like def convergence_rates(f, F, a, b, num_experiments=14): from math import log from numpy import zeros expected = F(b) - F(a) n = zeros(num_experiments, dtype=int) E = zeros(num_experiments) r = zeros(num_experiments-1) for i in range(num_experiments): n[i] = 2**(i+1) computed = trapezoidal(f, a, b, n[i]) E[i] = abs(expected - computed) if i > 0: r_im1 = -log(E[i]/E[i-1])/log(n[i]/n[i-1]) # Truncate to two decimals: r[i-1] = float(’{:.2f}’.format(r_im1)) return r Making a test function is a matter of choosingf,F,a, andb, and then checking thevalueof ri for the largest i: def test_trapezoidal_conv_rate(): """Check empirical convergence rates against the expected value 2.""" from math import exp v = lambda t: 3*(t**2)*exp(t**3) V = lambda t: exp(t**3) a = 1.1; b = 1.9 r = convergence_rates(v, V, a, b, 14) print(r) tol = 0.01 msg = str(r[-4:]) # show last 4 estimated rates assert (abs(r[-1]) - 2) < tol, msg Running the test shows that all ri, except the first one, equal the target limit 2 within twodecimals.Thisobservationsuggestsa toleranceof10−2. 6.7 DoubleandTripleIntegrals 6.7.1 TheMidpointRuleforaDoubleIntegral Givena double integralovera rectangulardomain [a,b]×[c,d], ∫ b a ∫ d c f(x,y)dydx,