Page - 64 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python

64 3 Computing Integrals # Integral by the trapezoidal method numerical = 0.5*v(a) + 0.5*v(b) for i in range(1, n): numerical += v(a + i*dt) numerical *= dt F = lambda t: exp(t**3) exact_value = F(b) - F(a) error = abs(exact_value - numerical) rel_error = (error/exact_value)*100 print ’n=%d: %.16f, error: %g’ % (n, numerical, error) Unfortunately, the twootherproblemsremainand theyare fundamental. Supposeyouwant to compute another integral, say R1:1 1 e x2dx. Howmuchdo weneed tochange in thepreviouscode tocompute thenewintegral?Not somuch: the formula forvmustbe replacedbyanewformula the limitsaandb the anti-derivativeV is not easily known2 andcanbeomitted, and thereforewe cannotwriteout theerror the notation should be changed to be alignedwith the newproblem, i.e.,t and dtchanged toxandh These changes are straightforward to implement, but they are scattered around in the program, a fact that requires us to be very careful sowedonot introduce new programmingerrorswhilewemodifythecode. It isalsoveryeasytoforgettomake a requiredchange. With the previous code intrapezoidal.py,we can compute the new integralR1:1 1 e x2dxwithout touching themathematicalalgorithm. Inan interactivesession (or inaprogram)wecan justdo >>> from trapezoidal import trapezoidal >>> from math import exp >>> trapezoidal(lambda x: exp(-x**2), -1, 1.1, 400) 1.5268823686123285 Whenyounowlookbackat the two solutions, theflat special-purposeprogram and the function-basedprogramwith the general-purpose functiontrapezoidal, youhopefullyrealize that implementingageneralmathematicalalgorithminagen- eral function requires somewhatmoreabstract thinking, but the resulting codecan beusedoverandoveragain. Essentially, ifyouapply theflat special-purposestyle, you have to retest the implementation of the algorithm after every change of the program. 2You cannot integrate e x2 by hand, but this particular integral is appearing so often in somany contexts that the integral is a special function, called theError function (http://en.wikipedia.org/ wiki/Error_function) andwritten erf.x/. In a code, you can call erf(x). The erf function is found in themathmodule.