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

74 3 Computing Integrals Note the importance of checking err against exactwith a tolerance: rounding errors from the arithmetics inside trapezoidalwill notmake the result exactly like thehand-computedone. Thesizeof the tolerance ishere set to10 14,which is a kindof all-roundvalue for computationswith numbersnot deviatingmuch from unity. Solvingaproblemwithoutnumerical errors Weknowthat the trapezoidal rule isexact for linear integrands.Choosingthe integral R4:4 1:2 .6x 4/dx as test case, the correspondingtest function for thisunit testmay look like def test_trapezoidal_linear(): """Check that linear functions are integrated exactly.""" f = lambda x: 6*x - 4 F = lambda x: 3*x**2 - 4*x # Anti-derivative a = 1.2; b = 4.4 expected = F(b) - F(a) tol = 1E-14 for n in 2, 20, 21: computed = trapezoidal(f, a, b, n) error = abs(expected - computed) success = error < tol msg = ’n=%d, err=%g’ % (n, error) assert success, msg Demonstrating correct convergence rates In the present examplewith integra- tion, it is known that the approximation errors in the trapezoidal rule are propor- tional ton 2,nbeing thenumberof subintervalsused in thecomposite rule. Computing convergence rates requires somewhat more tedious programming than the previous tests, but can be applied to more general integrands. The al- gorithmtypicallygoes like for i D0;1;2;:: :;q – ni D2iC1 – Compute integralwithni intervals – Compute theerrorEi – Estimate ri from(3.24) if i >0 Thecorrespondingcodemay look 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(float(n[i])/n[i-1]) r[i-1] = float(’%.2f’ % r_im1) # Truncate to two decimals return r