Principles of Numerical Calculations

85 Pages·2008·1.32 MB·English

dqbjvol1
2007/12/28
page1
Chapter 1
Principles of Numerical
Calculations
Itisalmostimpossibletoidentifyamathematical
theorynomatterhow“pure,”thathasnever
influencednumericalreasoning.
—B.J.C.BaxterandAriehIserles
1.1 Common Ideas and Concepts
Althoughnumericalmathematicshasbeenusedforcenturiesinoneformoranotherwithin
manyareasofscienceandindustry,1modernscientificcomputingusingelectroniccomput-
ershasitsorigininresearchanddevelopmentsduringtheSecondWorldWar. Inthelate
1940sandearly1950sthefoundationofnumericalanalysiswaslaidasaseparatediscipline
ofmathematics. Thenewcapabilityofperformingbillionsofarithmeticoperationscheaply
hasledtonewclassesofalgorithms,whichneedcarefulanalysistoensuretheiraccuracy
andstability.
Asarule,applicationsleadtomathematicalproblemswhichintheircompleteform
cannotbeconvenientlysolvedwithexactformulas,unlessonerestrictsoneselftospecial
casesorsimplifiedmodels. Inmanycases,onetherebyreducestheproblemtoasequenceof
linearproblems—forexample,linearsystemsofdifferentialequations. Suchanapproach
canquiteoftenleadtoconceptsandpointsofviewwhich,atleastqualitatively,canbeused
evenintheunreducedproblems.
Recentdevelopmentshaveenormouslyincreasedthescopeforusingnumericalmeth-
ods. Gainsinproblemsolvingcapabilitiesmeanthattodayonecantreatmuchmorecomplex
andlesssimplifiedproblemsthroughmassiveamountsofnumericalcalculations. Thishas
increasedtheinteractionofmathematicswithscienceandtechnology.
In most numerical methods, one applies a small number of general and relatively
simpleideas. Thesearethencombinedwithoneanotherinaninventivewayandwithsuch
1TheGreekmathematicianArchimedes(287–212B.C.),IsaacNewton(1642–1727),andCarlFriedrichGauss
(1777–1883)werepioneeringcontributorstonumericalmathematics.
1
Copyright ©2008 by the Society for Industrial and Applied Mathematics.
This electronic version is for personal use and may not be duplicated or distributed.

From "Numerical Methods in Scientific Computing, Volume 1" by Germund Dalquist and Åke Björck.
This book is available for purchase at www.siam.org/catalog.
dqbjvol1
2007/12/28
page2
2 Chapter1. PrinciplesofNumericalCalculations
knowledge of the given problem as one can obtain in other ways—for example, with the
methodsofmathematicalanalysis. Someknowledgeofthebackgroundoftheproblemis
alsoofvalue;amongotherthings,oneshouldtakeintoaccounttheordersofmagnitudeof
certainnumericaldataoftheproblem.
In this chapter we shall illustrate the use of some general ideas behind numerical
methods on some simple problems. These may occur as subproblems or computational
detailsoflargerproblems,thoughasaruletheymoreoftenoccurinalesspureformand
on a larger scale than they do here. When we present and analyze numerical methods,
tosomedegreeweusethesameapproachwhichwasfirstmentionedabove: westudyin
detail special cases and simplified situations, with the aim of uncovering more generally
applicableconceptsandpointsofviewwhichcanguideusinmoredifficultproblems.
Itisimportanttokeepinmindthatthesuccessofthemethodspresenteddependson
the smoothness properties of the functions involved. In this first survey we shall tacitly
assumethatthefunctionshaveasmanywell-behavedderivativesasareneeded.
1.1.1 Fixed-Point Iteration
One of the most frequently occurring ideas in numerical calculations is iteration (from
the Latin iterare, “to plow once again”) or successive approximation. Taken generally,
iterationmeanstherepetitionofapatternofactionorprocess. Iterationinthissenseoccurs,
forexample,intherepeatedapplicationofanumericalprocess—perhapsverycomplicated
anditselfcontainingmanyinstancesoftheuseofiterationinthesomewhatnarrowersense
tobedescribedbelow—inordertoimprovepreviousresults. Toillustrateamorespecific
useoftheideaofiteration,weconsidertheproblemofsolvinga(usually)nonlinearequation
oftheform
x =F(x), (1.1.1)
whereF isassumedtobeadifferentiablefunctionwhosevaluecanbecomputedforany
given value of a real variable x, within a certain interval. Using the method of iteration,
onestartswithaninitialapproximationx ,andcomputesthesequence
0
x =F(x ), x =F(x ), x =F(x ),.... (1.1.2)
1 0 2 1 3 2
Each computation of the type xn+1 = F(xn), n = 0,1,2,..., is called a fixed-point
iteration. Asngrows, wewouldlikethenumbersxn tobebetterandbetterestimatesof
thedesiredroot. Ifthesequence{xn}convergestoalimitingvalueα,thenwehave
α =nl→im∞xn+1 =nl→im∞F(xn)=F(α),
andthusx =αsatisfiestheequationx =F(x). Onecanthenstoptheiterationswhenthe
desiredaccuracyhasbeenattained.
AgeometricinterpretationoffixedpointiterationisshowninFigure1.1.1. Arootof
(1.1.1)isgivenbytheabscissa(andordinate)ofanintersectingpointofthecurvey =F(x)
andtheliney =x. Startingfromx ,thepointx =F(x )onthex-axisisobtainedbyfirst
0 1 0
drawing a horizontal line from the point (x ,F(x )) = (x ,x ) until it intersects the line
0 0 0 1
y = x inthepoint(x ,x ); fromtherewedrawaverticallineto (x ,F(x )) = (x ,x ),
1 1 1 1 1 2
Copyright ©2008 by the Society for Industrial and Applied Mathematics.
This electronic version is for personal use and may not be duplicated or distributed.

From "Numerical Methods in Scientific Computing, Volume 1" by Germund Dalquist and Åke Björck.
This book is available for purchase at www.siam.org/catalog.
dqbjvol1
2007/12/28
page3
1.1. CommonIdeasandConcepts 3
1 1
y = x y = x
0.8 0.8
y = F(x)
0.6 y = F(x) 0.6
0.4 0.4
0.2 0.2
x x x x x x x x
0 1 2 0 2 4 3 1
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
(a) 0 < F′(x) < 1 (b) −1 < F′(x) < 0
1 1
y = F(x) y = x
0.8 0.8
y = x
0.6 0.6
0.4 0.4
y = F(x)
0.2 0.2
x x x x x x x x x
3 2 1 0 4 2 0 1 3
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
(c) F′(x) > 1 (d) F′(x) < −1
Figure1.1.1. Geometricinterpretationofiterationxn+1 =F(xn).
and so on in a “staircase” pattern. In Figure1.1.1(a) it is obvious that the sequence {xn}
convergesmonotonicallytotherootα. Figure1.1.1(b)showsacasewhereF isadecreasing
function. Therewealsohaveconvergence,butnotmonotoneconvergence;thesuccessive
iterates xn lie alternately to the right and to the left of the root α. In this case the root is
bracketedbyanytwosuccessiveiterates.
Therearealsotwodivergentcases,exemplifiedbyFigures1.1.1(c)and(d). Onecan
see geometrically that the quantity, which determines the rate of convergence (or diver-
gence),istheslopeofthecurvey = F(x)intheneighborhoodoftheroot. Indeed,from
themeanvaluetheoremofcalculuswehave
x −α F(x )−F(α)
n+1 = n =F(cid:20)(ξ ),
x −α x −α n
n n
where ξn lies between xn and α. We see that if x0 is chosen sufficiently close to the root
(yetx (cid:30)= α),theiterationwillconvergeif|F(cid:20)(α)| < 1. Inthiscaseα iscalledapointof
0
attraction. Theconvergenceisfasterthesmaller|F(cid:20)(α)|is.
If|F(cid:20)(α)|>1,thenαisapointofrepulsionandtheiterationdiverges.
Copyright ©2008 by the Society for Industrial and Applied Mathematics.
This electronic version is for personal use and may not be duplicated or distributed.

From "Numerical Methods in Scientific Computing, Volume 1" by Germund Dalquist and Åke Björck.
This book is available for purchase at www.siam.org/catalog.
dqbjvol1
2007/12/28
page4
4 Chapter1. PrinciplesofNumericalCalculations
2.5
2
1.5
1
x x x
0.5 0 2 1
0.5 1 1.5 2 2.5
Figure1.1.2. Thefixed-pointiterationxn+1 =(xn+c/xn)/2,c=2,x0 =0.75.
Example1.1.1.
The square root of c > 0 satisfies the equation x2 = c, which can also be written
x =c/x orx =1/2(x+c/x). Thissuggeststhefixed-pointiteration
1
xn+1 = (xn+c/xn), n=1,2,..., (1.1.3)
2
whichisthewidelyusedHeron’srule.2 Thecurvey = F(x)isinthiscaseahyperbola
(seeFigure1.1.2).
From(1.1.3)follows
(cid:2) (cid:3) √
xn+1±√c= 21 xn±2√c+ xcn = (xn±2xn c)2;
thatis,
√ (cid:2) √ (cid:3)
x − c x − c 2
n+1 √ = n √ . (1.1.4)
x + c x + c
n+1 n
√
Wecantakeen = xxnn−+√cc tobeameas√ureoftheerr√orinxn. Then(1.1.4)readsen+1 = en2
anditfollowsthaten =e02n. If|x0− c|(cid:30)=|x0+ c|,thene0 <1andxn convergestoa
squarerootofcwhenn→∞. Notethattheiteration(1.1.3)canalsobeusedforcomplex
valuesofc.
Forc = 2andx = 1.5,wegetx = 1/2(1.5+2/1.5) = 15/12 = 1.4166666...,
0 1
and
x =1.414215686274, x =1.414213562375
2 3
√
(correct digits shown in boldface). This can be compared with the exact value 2 =
1.414213562373.... AscanbeseenfromFigure1.1.2,aroughvalueforx suffices. The
0
2Heronmadeimportantcontributionstogeometryandmechanics.HeisbelievedtohavelivedinAlexandria,
Egypt,duringthefirstcenturyA.D.
Copyright ©2008 by the Society for Industrial and Applied Mathematics.
This electronic version is for personal use and may not be duplicated or distributed.

From "Numerical Methods in Scientific Computing, Volume 1" by Germund Dalquist and Åke Björck.
This book is available for purchase at www.siam.org/catalog.
dqbjvol1
2007/12/28
page5
1.1. CommonIdeasandConcepts 5
√
rapidconvergenceisduetothefactthatforα = cwehave
F(cid:20)(α)=(1−c/α2)/2=0.
Onecaninfactshowthat √
|x − c|
lim n+1 √ =C
n→∞ |xn− c|2
for some constant 0 < C <∞, which is an example of what is known as quadratic
convergence. Roughly,ifxnhast correctdigits,thenxn+1willhaveatleast2t−1correct
digits.
The above iteration method is used quite generally on both pocket calculators and
largercomputersforcalculatingsquareroots.
Iterationisoneofthemostimportantaidsforpracticalaswellastheoreticaltreatment
ofbothlinearandnonlinearproblems. Oneverycommonapplicationofiterationistothe
solutionofsystemsofequations. Inthiscase{xn}isasequenceofvectors,andF isavector-
valuedfunction. Wheniterationisappliedtodifferentialequations,{xn}meansasequence
of functions, and F(x) means an expression in which integration or other operations on
functions may be involved. A number of other variations on the very general idea of
iterationwillbegiveninlaterchapters.
The form of (1.1.1) is frequently called the fixed-point form, since the root α is a
fixedpointofthemappingF. Anequationmaynotbegiveninthisformoriginally. One
has a certain amount of choice in the rewriting of an equation f(x) = 0 in fixed-point
form,andtherateofconvergencedependsverymuchonthischoice. Theequationx2 =c
canalsobewritten,forexample,asx = c/x. Theiterationformulaxn+1 = c/xn givesa
sequencewhichalternatesb√etweenx0(forevenn)andc/x0(foroddn)—thesequencedoes
notconvergeforanyx (cid:30)= c!
0
1.1.2 Newton’s Method
Letanequationbegivenintheformf(x)=0,andforanyk (cid:30)=0,set
F(x)=x+kf(x).
Then the equation x = F(x) is equivalent to the equation f(x) = 0. Since F(cid:20)(α) =
1+kf(cid:20)(α),weobtainthefastestconvergencefork =−1/f(cid:20)(α). Becauseαisnotknown,
thiscannotbeappliedliterally. Butifweusexnasanapproximation,thisleadstothechoice
F(x)=x−f(x)/f(cid:20)(x),ortheiteration
f(x )
xn+1 =xn− f(cid:20)(xn). (1.1.5)
n
ThisisthecelebratedNewton’smethod.3 Weshallderiveitinanotherwaybelow.
3IsaacNewton(1642–1727),Englishmathematician,astronomer,andphysicistinventedinfinitesimalcalculus
independentlyoftheGermanmathematicianandphilosopherGottfriedW.vonLeibniz(1646–1716).
Copyright ©2008 by the Society for Industrial and Applied Mathematics.
This electronic version is for personal use and may not be duplicated or distributed.

From "Numerical Methods in Scientific Computing, Volume 1" by Germund Dalquist and Åke Björck.
This book is available for purchase at www.siam.org/catalog.
dqbjvol1
2007/12/28
page6
6 Chapter1. PrinciplesofNumericalCalculations
Theequationx2 =ccanbewrittenintheformf(x)=x2−c=0. Newton’smethod
forthisequationbecomes
(cid:4) (cid:5)
x2−c 1 c
xn+1 =xn− n2xn = 2 xn+ xn , n=0,1,2,..., (1.1.6)
which is the fast method in Example 1.1.1. More generally, Newton’s method applied to
theequationf(x)=xp−c =0canbeusedtocomputec1/p,p =±1,±2,...,fromthe
iteration
xp−c
x =x − n .
n+1 n pxp−1
n
Thiscanbewrittenas
(cid:2) (cid:3)
xn+1 = p1 (p−1)xn+ xpc−1 = (−xnp)[(1−p)−cxn−p]. (1.1.7)
n
It is convenient to use the first expression in (1.1.7) when p > 0 and the second when
p√<√0. Withp√= 2,3,and−2,respectively,thisiterationformulaisusedforcalculating
c, 3c,and1/ c. Also1/c,(p =−1)canbecomputedbytheiteration
xn+1 =xn+xn(1−cxn)=xn(2−cxn),
usingonlymultiplicationandaddition. Insomeearlycomputers,whichlackeddivisionin
hardware,thisiterationwasusedtoimplementdivision,i.e.,b/cwascomputedasb(1/c).
Example1.1.2.
WewanttoconstructanalgorithmbasedonNewton’smethodfortheefficientcalcu-
lationofthesquarerootofanygivenfloating-pointnumbera. Ifwefirstshiftthemantissa
sothattheexponentbecomeseven,a =c·22e,and1/2≤c<2;then
√ √
a = c·2e.
W√eneedonlyconsiderthereducedrange1/2≤c≤1sincefor1<c≤2wecancompute
1/candinvert.4
Tofindaninitialapproximationx√0 tostarttheNewtoniterationswhen1/2≤c <1,
wecanuselinearinterpolationofx = cbetweentheendpoints1/2,1,giving
√
x (c)= 2(1−c)+2(c−1/2)
0
√
( 2isprecomputed√). Theiterationthenproceedsusing(1.1.6). √
Forc = 3/4( c = 0.86602540378444)theresultisx = ( 2+2)/4and(correct
0
digitsareinboldface)
x =0.85355339059327, x =0.86611652351682,
0 1
x =0.86602540857756, x =0.86602540378444,
2 3
4Sincedivisionisusuallymuchslowerthanadditionandmultiplication,thismaynotbeoptimal.
Copyright ©2008 by the Society for Industrial and Applied Mathematics.
This electronic version is for personal use and may not be duplicated or distributed.

From "Numerical Methods in Scientific Computing, Volume 1" by Germund Dalquist and Åke Björck.
This book is available for purchase at www.siam.org/catalog.
dqbjvol1
2007/12/28
page7
1.1. CommonIdeasandConcepts 7
The quadratic rate of convergence is apparent. Three iterations suffice to give about 16
digitsofaccuracyforallx ∈[1/2,1].
Newton’smethodisbasedonlinearization. Thismeansthatlocally,i.e.,inasmall
neighborhoodofapoint,amorecomplicatedfunctionisapproximatedwithalinearfunc-
tion. Inthesolutionoftheequationf(x) = 0,thismeansgeometricallythatweseekthe
intersection point between the x-axis and the curve y = f(x); see Figure1.1.3. Assume
x
0
x x
2 1
Figure1.1.3. GeometricinterpretationofNewton’smethod.
thatwehaveanapproximatingvaluex totheroot. Wethenapproximatethecurvewithits
0
tangentatthepoint(x ,f(x )). Letx betheabscissaofthepointofintersectionbetween
0 0 1
thex-axisandthetangent. Sincetheequationforthetangentreads
y−f(x )=f(cid:20)(x )(x−x ),
0 0 0
bysettingy =0weobtaintheapproximation
x =x −f(x )/f(cid:20)(x ).
1 0 0 0
In many cases x will have about twice as many correct digits as x . But if x is a poor
1 0 0
approximationandf(x)farfromlinear,thenitispossiblethatx willbeaworseapproxi-
1
mationthanx .
0
Ifwecombinetheideasofiterationandlinearization,thatis,substitutexnforx0and
xn+1 forx1,werediscoverNewton’smethodmentionedearlier. Ifx0 iscloseenoughtoα,
theiterationswillconvergerapidly(seeFigure1.1.3),buttherearealsocasesofdivergence.
An alternative to drawing the tangent to approximate a curve locally with a linear
function is to choose two neighboring points on the curve and to approximate the curve
with the secant which joins the two points; see Figure1.1.4. The secant method for the
solutionofnonlinearequationsisbasedonthisapproximation. Thismethod,whichpreceded
Newton’smethod,isdiscussedinmoredetailinSec.6.3.1.
Newton’smethodcanbegeneralizedtoyieldaquadraticallyconvergentmethodfor
solvingasystemofnonlinearequations
fi(x1,x2,...,xn)=0, i =1:n. (1.1.8)
Copyright ©2008 by the Society for Industrial and Applied Mathematics.
This electronic version is for personal use and may not be duplicated or distributed.

From "Numerical Methods in Scientific Computing, Volume 1" by Germund Dalquist and Åke Björck.
This book is available for purchase at www.siam.org/catalog.
dqbjvol1
2007/12/28
page8
8 Chapter1. PrinciplesofNumericalCalculations
x x x
0 2 3
x
1
Figure1.1.4. Geometricinterpretationofthesecantmethod.
Suchsystemsariseinmanydifferentcontextsinscientificcomputing. Importantexamples
are the solution of systems of differential equations and optimization problems. We can
write(1.1.8)asf(x)=0,wheref andxarevectorsinRn. Thevector-valuedfunctionf is
saidtobedifferentiableatthepointxifeachcomponentisdifferentiablewithrespecttoall
thevariables. Thematrixofpartialderivativesoff withrespecttox,
 
∂f1 ... ∂f1
∂x ∂x
J(x)=f(cid:20)(x)= ...1 ...n ∈Rn×n, (1.1.9)
∂fn ... ∂fn
∂x ∂x
1 n
iscalledtheJacobianoff.
Letxkbethecurrentapproximatesolutionandassumethatthematrixf(cid:20)(xk)isnonsin-
gular. TheninNewton’smethodforthesystem(1.1.8),thenextiteratexn+1isdetermined
fromtheuniquesolutiontothesystemoflinearequations
J(xk)(xn+1−xk)=−f(xk). (1.1.10)
Thelinearsystem(1.1.10)canbesolvedbycomputingtheLUfactorizationofthematrix
J(x);seeSec.1.3.2.
EachstepofNewton’smethodrequirestheevaluationofthen2entriesoftheJacobian
matrix J(xk). This may be a time consuming task if n is large. If either the iterates or
the Jacobian matrix are not changing too rapidly, it is possible to reevaluate J(xk) only
occasionallyandusethesameJacobianinseveralsteps. Thishasthefurtheradvantagethat
oncewehavecomputedtheLUfactorizationoftheJacobianmatrix,thelinearsystemcan
besolvedinonlyO(n2)arithmeticoperations;seeSec.1.3.2.
Example1.1.3.
ThefollowingexampleillustratesthequadraticconvergenceofNewton’smethodfor
simpleroots. Thenonlinearsystem
x2+y2−4x =0,
y2+2x−2=0
Copyright ©2008 by the Society for Industrial and Applied Mathematics.
This electronic version is for personal use and may not be duplicated or distributed.

From "Numerical Methods in Scientific Computing, Volume 1" by Germund Dalquist and Åke Björck.
This book is available for purchase at www.siam.org/catalog.
dqbjvol1
2007/12/28
page9
1.1. CommonIdeasandConcepts 9
hasasolutionclosetox =0.5,y =1. TheJacobianmatrixis
0 0
(cid:2) (cid:3)
2x−4 2y
J(x,y)= ,
2 2y
andNewton’smethodbecomes
(cid:2) (cid:3) (cid:2) (cid:3) (cid:2) (cid:3)
xk+1 = xk −J(x ,y )−1 xk2+yk2−4xn .
yk+1 yk k k yk2+2xk−2
Wegetthefollowingresults:
k x y
k k
1 0.35 1.15
2 0.35424528301887 1.13652584085316
3 0.35424868893322 1.13644297217273
4 0.35424868893541 1.13644296914943
Alldigitsarecorrectinthelastiteration. Thequadraticconvergenceisobvious;thenumber
ofcorrectdigitsapproximatelydoublesineachiteration.
Often,themaindifficultyinsolvinganonlinearsystemistofindasufficientlygood
starting point for the Newton iterations. Techniques for modifying Newton’s method to
ensure global convergence are therefore important in several dimensions. These must
include techniques for coping with ill-conditioned or even singular Jacobian matrices at
intermediatepoints. SuchtechniqueswillbediscussedinVolumeII.
1.1.3 Linearization and Extrapolation
Thesecantapproximationisusefulinmanyothercontexts;forinstance,itisgenerallyused
whenone“readsbetweenthelines”orinterpolatesinatableofnumericalvalues. Inthiscase
thesecantapproximationiscalledlinearinterpolation. Whenthesecantapproximationis
usedinnumericalintegration,i.e.,intheapproximatecalculationofadefiniteintegral,
(cid:12)
b
I = y(x)dx, (1.1.11)
a
(seeFigure1.1.5)itiscalledthetrapezoidalrule. Withthismethod,theareabetweenthe
curvey = y(x)andthex-axisisapproximatedwiththesumT(h)oftheareasofaseries
ofparalleltrapezoids. UsingthenotationofFigure1.1.5,wehave
1(cid:13)n−1 b−a
T(h)=h2 (yi +yi+1), h= n . (1.1.12)
i=0
(In the figure, n = 4.) We shall show in a later chapter that the error is very nearly
proportional to h2 when h is small. One can then, in principle, attain arbitrarily high
accuracybychoosinghsufficientlysmall. Butthecomputationalworkinvolvedisroughly
Copyright ©2008 by the Society for Industrial and Applied Mathematics.
This electronic version is for personal use and may not be duplicated or distributed.

From "Numerical Methods in Scientific Computing, Volume 1" by Germund Dalquist and Åke Björck.
This book is available for purchase at www.siam.org/catalog.
dqbjvol1
2007/12/28
page10
10 Chapter1. PrinciplesofNumericalCalculations
y
0
y
4
y
1
y2 y3
a b
Figure1.1.5. Numericalintegrationbythetrapezoidalrule(n=4).
proportional to the number of points where y(x) must be computed, and thus inversely
proportional to h. Hence the computational work grows rapidly as one demands higher
accuracy(smallerh).
Numericalintegrationisafairlycommonproblembecauseonlyseldomcanthe“prim-
itive”functionbeanalyticallycalculatedinafiniteexpressioncontainingonlyelementary
functions. Itisnotpossibleforsuchsimplefunctionsasex2 or(sinx)/x. Inordertoobtain
higheraccuracywithsignificantlylessworkthanthetrapezoidalrulerequires,onecanuse
oneofthefollowingtwoimportantideas:
(a) localapproximationoftheintegrandwithapolynomialofhigherdegree,orwitha
functionofsomeotherclass,forwhichoneknowstheprimitivefunction;
(b) computationwiththetrapezoidalruleforseveralvaluesofhandthenextrapolation
to h = 0, the so-called Richardson extrapolation5 or deferred approach to the
limit,withtheuseofgeneralresultsconcerningthedependenceoftheerroronh.
Thetechnicaldetailsforthevariouswaysofapproximatingafunctionwithapoly-
nomial, including Taylor expansions, interpolation, and the method of least squares, are
treatedinlaterchapters.
Theextrapolationtothelimitcaneasilybeappliedtonumericalintegrationwiththe
trapezoidalrule. Aswasmentionedpreviously, thetrapezoidalapproximation(1.1.12)to
the integral has an error approximately proportional to the square of the step size. Thus,
usingtwostepsizes,hand2h,onehas
T(h)−I ≈kh2, T(2h)−I ≈k(2h)2,
andhence4(T(h)−I)≈T(2h)−I,fromwhichitfollowsthat
1 1
I ≈ (4T(h)−T(2h))=T(h)+ (T(h)−T(2h)).
3 3
5LewisFryRichardson(1881–1953)studiedmathematics,physics,chemistry,botany,andzoology. Hegrad-
uatedfromKing’sCollege,Cambridgein1903. Hewasthefirst(1922)toattempttoapplythemethodoffinite
differencestoweatherprediction,longbeforethecomputerage!
Copyright ©2008 by the Society for Industrial and Applied Mathematics.
This electronic version is for personal use and may not be duplicated or distributed.

From "Numerical Methods in Scientific Computing, Volume 1" by Germund Dalquist and Åke Björck.
This book is available for purchase at www.siam.org/catalog.