Finite Di erence Methods for Di erential Equations

1 FiniteDi erenceMethodsforDi erentialEquationsRandallJ. LeVequeDRAFTVERSION forusein thecourseAMath585{586 University of WashingtonVersionof September,2005 WARNING:Thesenotesareincompleteandmay in my R. J. LeVeque,1998{20052c LeVeque,2004|University of Washington|AMath585{6 NotesContentsIBasicText11 Finitedi .. nitedi erenceapproximations..82 .. simple nitedi erencemethod .. in the2-norm.. 'sfunctionsandmax-normstability .. generallinearsecondorderequation.. thenonlinearBVP.. re nement .. erencing..37iiiCONTENTS3 .. ve-point stencilfortheLaplacian.. Laplacian..484 .. basedonpolynomialinterpolation.. niteelement method .. spacedimensions..605 Iterative Methods .. matrixsplittingmethods.. convergence.. Methods andconjugategradients.. of steepestdescent .. direction.. algorithm.. CG..866 TheInitialValueProblemforODE' .. solutions.. Equations .. canceof theLipschitzconstant .. LeVeque|AMath585{6 Notesiii7 Zero-Stability.}}}}

2 'sprinciple.. 'smethod onlinearproblems.. stability forBVP's.. 'smethod onnonlinearproblems.. of linearmultistepmethods.. erenceequations.. 1188 AbsoluteStability zero-stablemethod .. regionsforLMMs.. as one-stepmethods ona system.. ordinarydi erentialequations.. stepsize.. 1349 Sti culties.. sti ness.. forsti problems.. 14110 cationof di erentialequations.. PDEsfromconservationprinciples.. usion.. usionequations.. 148ivCONTENTS11 FourierAnalysisof .. di erentialequations.. waves.. 15212 Di accuracy.. of Linesdiscretizations.. theory.. nessof theheatequation.. theory.. 16813 .. method .. analysis.. analysis.. edEquations.. waves.. packets.. forhyperbolicsystems.. 19014 erences.. LeVeque|AMath585{6 Notesv15 .. fractionalstepmethods.. 207 IIAppendicesA{1A1 MeasuringErrorsA{ a scalarvalue.. A{ .. A{ error.. A{ \Big-oh"and\little-oh"notation.. A{ vectors.}}}}}}}

3 A{ .. A{ .. A{ functions.. A{ gridfunctions.. A{ .. A{7A2 Estimatingerrorsin numericalsolutionsA{ .. A{ ne-gridsolution.. A{ .. A{11A3 EigenvaluesandinnerproductnormsA{ transformations.. A{ .. A{ .. A{ .. A{ .. A{ .. A{ matrices.. A{ .. A{ .. A{ .. A{ .. A{25A4 MatrixpowersandexponentialsA{ matrics.. A{ .. A{ .. A{ non-normality.. A{ .. A{ matricesandtheKreissMatrixTheorem.. A{35A5 LinearDi erentialandDi erenceEquationsA{ erentialequations.. A{ erenceequations.. A{ .. A{40viCONTENTSc LeVeque,2004|University of Washington|AMath585{6 NotesPartIBasicText1c LeVeque,2004|University of Washington|AMath585{6 NotesChapter1 Finitedi erenceapproximationsOurgoalis to approximatesolutionsto di erentialequations, , to nda function(orsomediscreteapproximationto thisfunction)which satis esa given relationshipbetweenvariousof itsderivatives onsomegivenregionof spaceand/ortime,alongwithsomeboundarycon ditionsalongtheedgesof generalthisis a di cultproblemandonlyrarelycanananalyticfor mulabe nitedi erencemethod proceedsby replacingthederivatives in thedi erentialequationsby nitedi a largealgebraicsystemof equationsto be solvedinplaceof thedi erentialequation,somethingthatis easilysolvedona ,we rstconsiderthemorebasicquestionof how we canapproximatethederivatives of a knownfunctionby nitedi erenceformulasbasedonlyonvaluesof thefunctionitselfat basisforthelaterdevelopment of nitedi erencemethodsforsolvingdi erentialequations,thisallowsus to investigateseveralkey conceptssuch as theorderofaccuracyof anapproximationin (x)}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}

4 Represent a functionof onevariablethat,unlessotherwisestated,wi llalways be assumedto be smooth,meaningthatwe candi erentiatethefunctionseveraltimesandeach derivative is awell-de nedboundedfunctionover aninterval containinga particularpoint of interest want to approximateu0( x) by a nitedi erenceapproximationbasedonlyonvaluesofua t a nitenumber of points near x. Oneobviouschoicewouldbe to useD+u( x) u( x+h) u( x)h( )forsomesmallvalueofh. Thisis motivatedby thestandardde nitionof thederivative as thelimitingvalueof thisexpressionash! +u( x) is theslope of thelineinterpolatinguat thepoints xand x+h( ).Theexpression( )is aone-sidedapproximationtou0sinceuis evaluatedonlyat valuesofx u( x) u( x) u( x h)h:( )Each of theseformulasgives a rstorderaccurateapproximationtou0( x), meaningthatthesizeof theerroris is to usethecentered approximationD0u( x) u( x+h) u( x h)2h=12(D+u( x) +D u( x)):( )Thisis theslope of thelineinterpolatinguat x hand x+h, andis simplytheaverageof thetwoone-sidedapproximationsde shouldbe clearthatwe wouldexpectD0u( x) to give a betterapproximationthaneitherof factthisgives a34 Finitedi erenceapproximationsPSfragreplacements x h x x+hu(x)slopeu0( x)slopeD+u( x)slopeD u( x)slopeD0u( x) :Variousapproximationstou0( x) interpretedas theslope of secant :Errorsin various nitedi erenceapproximationstou0( x).

5 HD+ | theerroris proportionaltoh2andhenceis much smallerthantheerrorin a rstorderapproximationwhenhis ,forexampleD3u( x) 16h[2u( x+h) + 3u( x) 6u( x h) +u( x 2h)]:( )It may notbe clearwherethiscamefromor why it shouldapproximateu0at all,butin factit turnsoutto be a thirdorderaccurateapproximation|theerror is proportionaltoh3whenhis rstgoalis to developsystematicways to derive such formulasandto analyzetheiraccuracyandrelative willlook at a typicalexampleof how theerrorsin (x) = sin(x) and x= 1, so we aretryingto approximateu0(1)= cos(1)=0 ( x) u0( x) forvariousvaluesofhforeach of seethatD+uandD ubehave similarlythoughoneexhibitsan errorthatis roughlythenegativeof , theaverageof thetwo, hasanerrorthatis much seethatD+u( x) u0( x) 0:42hD0u( x) u0( x) 0:09h2D3u( x) u0( x) 0 LeVeque|AMath585{6 Notes510 310 210 110 1010 810 610 410 2 PSfragreplacementsD+ :TheerrorsinDu( x) rmingthatthesemethods are rstorder,secondorder,andthirdorder, a good way to ploterrorswhenwe expectthemto behave like somepower ofh, sinceif theerrorE(h) behaves likeE(h) ChpthenlogjE(h)j logjCj+plogh:Soona log-logscaletheerrorbehaves linearlywitha slope thatis equaltop, theorderof to analyzingtheerrorin a nitedi erenceapproximationis to expandeach ofthefunctionvaluesofuin aTaylorseriesaboutthepoint x, ,u( x+h)=u( x) +hu0( x) +12h2u00( x) +16h3u000( x) +O(h4)( )u( x h)=u( x) hu0( x) +12h2u00( x) 16h3u000( x) +O(h4)( )Theseexpansionsarevalidprovidedthatuis su \big-oh"notationO(h4) areadvisedto AppendixA1at thispoint sincethisnotationwillbe heavilyusedanda proper understandingof itsuseis ( )allowsus to computethatD+u( x) =u( x+h) u( x)h=u0( x) +12hu00( x) +16h2u000( x) +O(h3):Recallthat xis a xedpoint so thatu00( x); u000( x), etc.}

6 ,are xedconstants independent ofh. Theydependonuof course,butthefunctionis also xedas we cientlysmall,theerrorwillbe dominatedby the rstterm12hu00( x) andalltheothertermswillbe negligiblecomparedto thisterm,so we expecttheerrorto behave roughlylike a constanttimesh, wheretheconstant hasthevalue12u00( x).Notethatin ,whereu(x) = sinx, we have12u00(1)= 0:4207355which agreeswiththebehaviorseenin erenceapproximationsSimilarly, from( )we cancomputethattheerrorinD u( x) isD u( x) u0( x) = 12hu00( x) +16h2u000( x) +O(h3)which ( )and( )showsthatu( x+h) u( x h) = 2hu0( x) +13h3u000( x) +O(h5)so thatD0u( x) u0( x) =16h2u000( x) +O(h4):( )Thiscon rmsthesecondorderaccuracyof thisapproximationandagainagreeswithwhati s ,sincein thecontextof have16u000( x) = 16cos(1)= 0:09005038:Notethatallof theoddordertermsdropoutof theTaylorseriesexpansion( )forD0u( x). Thisistypicalwithcenteredapproximationsa ndtypicallyleadsto a orderto analyzeD3uwe needto alsoexpandu( x 2h) asu( x 2h) =u( x) 2hu0( x) +12(2h)2u00( x) 16(2h)3u000( x) +O(h4):( )Combiningthiswith( )and( )showsthatD3u( x) =u0( x) +112h3u0000( x) +O(h4):( ) nitedi erenceapproximationsSupposewe want to derive a nitedi erenceapproximationtou0( x) basedonsomegiven setof canuseTaylorseriesto derive an appropriateformula,usingthemethod of undetermined coe want a one-sidedapproximationtou0( x) basedonu( x); u( x h) andu( x 2h), of theformD2u( x) =au( x) +bu( x h) +cu( x 2h):( )We candeterminethecoe cientsa; b, andcto give thebestpossibleaccuracyby expandingin ( )and( )in ( )givesD2u( x)=(a+b+c)u( x) (b+ 2c)hu0( x) +12(b+ 4c)h2u00( x) 16(b+ 8c)h3u000( x) + :If thisis goingto agreewithu0( x) to highorderthenwe needa+b+c=0b+ 2c= 1=h( )b+ 4c=0We might like to requirethathigherordercoe cients be zeroas well,butsincethereareonlythreeunknownsa; b.

7 Andcwe cannotin generalhope to satisfymorethanthreesuch ( )givesa= 3=2h b= 2=hc= 1= LeVeque|AMath585{6 Notes7so thattheformulaisD2u( x) =12h[3u( x) 4u( x h) +u( x 2h)]:( )Theerrorin thisapproximationis clearlyD2u( x) u0( x)= 16(b+ 8c)h3u000( x) + =112h2u000( x) +O(h3) to derive thesame nitedi is to approximatethefunctionu(x) by somepolynomialp(x) andthenusep0( x) as anapproximationtou0( x).If wedeterminethepolynomialby interpolatinguat anappropriatesetof points,thenwe obtainthesame nitedi erencemethods as derive themethod of thisway, letp(x) be thequadraticpolynomialthatinterpolatesua t x, x hand x 2handthencomputep0( x). Theresultis exactly( ). thesecondderivativeu00(x) canbe obtainedin givenbyD2u( x)=1h2[u( x h) 2u( x) +u( x+h)]=u00( x) +12h2u0000( x) +O(h4):Again,sincethisis a symmetriccenteredapproximationallof obtainedby themethod of undeterminedcoe cients,or alternativelybycomputingthesecondderivat ive of thequadraticpolynomialinterpolatingu(x) at x h; xand x+ to derive approximationsto higherorderderivatives is by repeatedlyapplying rstorderdi thesecondderivative is thederivative ofu0, we canviewD2u( x) as beingadi erenceof rstdi fact,D2u( x) =D+D u( x)sinceD+(D u( x))=1h[D u( x+h) D u( x)]=1h u( x+h) u( x)h u( x) u( x h)h =D2u( x):Alternatively,D2( x) =D D+u( x) or we canalsoviewit as a centereddi erenceof centereddi erences,if we usea stepsizeh=2 in each centeredapproximationto the rstderivative.}

8 If we de ne^D0u(x) =1h(u(x+h=2) u(x h=2))thenwe ndthat^D0(^D0u( x)) =1h u( x+h) u( x)h u( x) u( x h)h =D2u( x):8 Finitedi erenceapproximationstohigherorderderivat ives canalsobe obtainedusingany of erencingapproximationstolower orderderivatives is ,herearetwo di erent approximationstou000( x). The rstoneisuncenteredand rstorderaccurate:D+D2u( x)=1h3(u( x+ 2h) 3u( x+h) + 3u( x) u( x h))=u000( x) +12hu0000( x) +O(h2):Thenextapproximationis centeredandsecondorderaccurate:D0D+D u( x)=12h3(u( x+ 2h) 2u( x+h) + 2u( x h) u( x 2h))=u000( x) +14h2u00000( x) +O(h4):Finitedi erenceapproximationsof thesortderivedabove arethebasicbuildingblocksof nitedi erencemethods forsolvingdi :PSfragreplacementsh1h2h3x1x2x3x41. Usepolynomialinterpolationto derivea nitedi erence approximationforu00(x2)thatis asaccurateas possibleforsmoothfunctionsu, based on thefourvaluesU1=u(x1); :::,U4=u(x4).Giveanexpressionforthedomin anttermin Verifyyourexpressionfortheerror by testingyourformulawitha speci cfunctionandvariousvaluesofh1; h2; Canyoude nean\orderof accuracy"foryourmethod in termsofh= max(h1; h2; h3)?

9 Togeta betterfeel forhowtheerror behavesas thegridgets ner, largenumber (say500)of di erentvaluesofHspanningtwoor three ordersof magnitude,chooseh1; h2, andh3as randomnumbersin theinterval[0;H]andcomputetheerror in a log-log plotto geta scatterplotof thebehaviorasH!0. (Note:inmatlabthecommandh = H * rand(1)will produce a singlerandomnumberuniformlydistributed in therange[0,H].) Of coursetheseerrorswill notlie exactlyon a straightlinesince thevaluesofhkmayvaryquitea lot evenforH's thatare nearby,butyoumightexpecttheupper limitof theerror to Estimatethe\orderof accuracy"by doinga leastsquares t of theformlog(E(H)) =K+plog(H) LeVeque|AMath585{6 Notes9to determineKandpbased thatthiscanbe doneby solvingthefollowinglinear systemin theleastsquares sense:266641log(H1)1log(H2)..1log(H500)3 7775 Kp =26664log(E(H1))log(E(H2))..log(E(H500)) 37775:Inmatlaba rectangularsystemAx=bcanbe solved in theleastsquares sensebyx = A\ of undetermined coe cientsto nda fourth-orderaccurate nitedi er-ence approximationtou00(x)based on5 equally spaced points,u00(x) =c 2u(x 2h) +c 1u(x h) +c0u(x) +c1u(x+h) +c2u(x+ 2h) +O(h4):Testyourformulaonsomesmoothfuncti onto verifythatit givestheexpected erenceapproximationsc LeVeque,2004|University of Washington|AMath585{6 NotesChapter2 BoundaryValueProblemsWe will rstconsiderordinarydi erentialequationsthatareposedonsomeinter vala < x < b,togetherwithsomeboundaryconditionsatea ch endof willextendthistomorethanonespacedimensio n,andstudyellipticpartialdi erentialequationsthatareposedin someregionof theplaneorthree-dimensionalspace,andares olvedsubjecttosomeboundaryconditionsspec ifyingthesolutionand/oritsderivatives aroundtheboundaryof thesetwo chaptersaregenerallysteadystateproblemsi n which thesolutionvariesonlywiththespatialcoord inatesbutnotwithtime.}}

10 ( casewhere[a;b]is a timeinterval ratherthananinterval in space.)Steady-stateproblemsareoftenassoc iatedwithsometime-dependent problemthatdescribes thedynamicbehavior,andthe2-point boundaryvalueproblemor ellipticequationresultsfromconsideringth especialcasewherethesolutionis steadyin time,andhencethetime-derivative termsareequaltozero, speci cexample,considerthe ow of heatin a rod madeoutof someheat-conductingmaterial,subjectto someexternalheatsourcealongitslengthands omeboundaryconditionsat each assumethatthematerialproperties,theiniti altemperaturedistribution,andthesourceva ryonlywithx, thedistancealongthelength,andnotacrossan y cross-section,thenwe expectthetemperaturedistributionat any timeto varyonlywithxandwe canmodelthiswitha di erentialequationin varywithtime,we letu(x;t) denotethetemperatureatpointxat timet, wherea < x < balongsome nitelengthof thengovernedby theheat equationut(x;t) = ( (x)ux(x;t))x+ (x;t)( )where (x) is thecoe cient of heatconduction,which may varywithx, and (x;t) is theheatsource(orsink,if <0).