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

70 3 Computing Integrals Solving a problemwithout numerical errors The best unit tests for numerical algorithmsinvolvemathematicalproblemswhereweknowthenumerical resultbe- forehand. Usually, numerical results contain unknown approximation errors, so knowing the numerical result implies that we have a problemwhere the approx- imation errors vanish. This feature may be present in very simplemathematical problems. For example, the trapezoidal method is exact for integration of linear functionsf.x/ D axCb. We can therefore pick some linear function and con- struct a test function that checks equality between the exact analytical expression for the integral and thenumbercomputedby the implementationof the trapezoidal method. Aspecific test casecanbe R4:4 1:2 .6x 4/dx. This integral involvesan“arbitrary” interval Œ1:2;4:4 andan“arbitrary” linear functionf.x/D6x 4. By“arbitrary” wemeanexpressionswhereweavoid the specialnumbers0and1since thesehave specialproperties inarithmeticoperations (e.g., forgetting tomultiply is equivalent tomultiplyingby1,and forgetting toadd is equivalent toadding0). Demonstratingcorrect convergencerates Normally, unit testsmustbebasedon problemswhere thenumerical approximationerrors inour implementation remain unknown. However,we often knowormay assume a certain asymptotic behavior of the error. Wecando someexperimental runswith the test problem R1 0 3t2et 3 dt wheren is doubled in each run: n D 4;8;16. The corresponding errors are then 12%,3%and0.77%,respectively.Thesenumbersindicate that theerrorisroughly reducedbya factorof4whendoublingn. Thus, theerrorconverges tozeroasn 2 andwe say that the convergence rate is 2. In fact, this result can also be shown mathematically for the trapezoidal andmidpoint method. Numerical integration methodsusuallyhaveanerror thatconvergetozeroasn p forsomep thatdepends on the method. With such a result, it does not matter if we do not knowwhat the actual approximation error is: we knowatwhat rate it is reduced, so running the implementation for two ormore differentn valueswill put us in a position to measure theexpected rateandsee if it is achieved. The idea of a corresponding unit test is then to run the algorithm for some n values, compute the error (the absolute value of the difference between the exact analytical resultandtheoneproducedbythenumericalmethod),andcheckthat the error has approximately correct asymptotic behavior, i.e., that the error is propor- tional ton 2 in caseof the trapezoidalandmidpointmethod. Let us developamoreprecisemethod for such unit tests basedonconvergence rates.Weassume that theerrorE dependsonnaccording to E DCnr; whereC is an unknown constant and r is the convergence rate. Consider a set of experimentswith variousn: n0;n1;n2;:: :;nq. We compute the corresponding errorsE0;:::;Eq. For twoconsecutive experiments, number i and i 1,wehave theerrormodel Ei DCnri; (3.22) Ei 1 DCnri 1 : (3.23)