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

6.2 TheCompositeTrapezoidalRule 141 v_sum = 0 for i in range(1, n, 1): t = a + i*dt v_sum = v_sum + v(t) numerical = dt*(0.5*v(a) + v_sum + 0.5*v(b)) V = lambda t: exp(t**3) exact_value = V(b) - V(a) error = abs(exact_value - numerical) rel_error = (error/exact_value)*100 print(’n={:d}: {:.16f}, error: {:g}’.format(n, numerical, error)) Unfortunately, the two other problems (1. and 3.) remain, and they are funda- mental. ComputingAnotherIntegral Supposeyounextwanttocomputeanotherintegral, say ∫1.1 −1 e −x2dx, using the previous specific implementation as your “starting point”.Whatchangesare required in thecode then? Firstofall, ananti-derivativecannot (easily)be found2,3 for thisnewintegrand, so we drop computing the integration error, and must remove the corresponding code lines. Inaddition, • thenotationshouldbechangedtofit thenewproblem.Thus,tanddt shouldbe replacedbyxandh.Also, the integrandis (most likely)notavelocityanymore, so thenamev shouldbechangedwith,e.g.,f.Similarly,v_sumshouldratherbe f_sum then. • the formula forv (orf)mustbe replacedbya newformula • the limitsaandbmust bechanged These changes are straightforward to implement, but they are scattered around in the program, a fact that requires us to be very careful so we do not introduce new programmingerrorswhile we modify the code. It is also veryeasy to forgetone or twoof the requiredchanges. For thesakeofcomparison,wemightseehoweasy it is to ratheruseourgeneral implementation in trapezoidal.py for the task. With the following interactive session, it should be clear that this implementation allows us to compute the new integral ∫1.1 −1 e −x2dxwithouttouchingtheimplementedmathematicalalgorithm!We cansimplydo: In [1]: from trapezoidal import trapezoidal # ...general implementation In [2]: from math import exp In [3]: trapezoidal(lambda x: exp(-x**2), -1, 1.1, 400) Out[3]: 1.5268823686123285 2 You cannot integrate e−x2 by hand, but this particular integral is appearing so often in so many contexts that the integral is a special function, called the Error function and written erf(x). In a code, you can callerf(x).Theerf function is found in themathmodule. 3 http://en.wikipedia.org/wiki/Error_function.