Computational physics

48 Pages·2002·0.316 MB·English

by  MacKinnon A.

Computational Physics — 3rd/4th Year Option
AngusMacKinnon
September26,2002
Computational Physics
Course Description
Theuseofcomputersinphysics,aswellasmostotherbranchesofscienceandengineering,hasincreased
manytimesalongwiththerapiddevelopmentoffasterandcheaperhardware.Thiscourseaimstogivethe
studentathoroughgroundinginthemaincomputationaltechniquesusedinmodernphysics. Itisparticu-
larlyimportantinthiscoursethatthestudentsshouldlearnbydoing.Thecourseisthereforedesignedsuch
that a significant fraction ofthe students’time isspent actually programmingspecific physicalproblems
ratherthanlearningabstracttechniques.
T(cid:0) hecoursewillcoverproblemsin4(5)broadsections:
(cid:0) Ordinarydifferentialequations,suchasthoseofclassicalmechanics.
Partialdifferentialequations,suchasMaxwell’sequationsandtheDiffusionandSchro¨dingerequa-
(cid:0) tions.
Matrixmethods,suchassystemsofequationsandeigenvalueproblemsappliedtoPoisson’sequation
(cid:0) andelectronicstructurecalculations.
MonteCarloandothersimulationmethods,suchastheMetropolisalgorithmandmoleculardynam-
(cid:0) ics.
(Iftimepermits:)ComputerAlgebra;anintroductionusingMapletotheusesandabusesofalgebraic
computinginphysics.
Thisisnotashortcourseincomputingscience,norinprogramming.Itfocusesspecificallyonmethods
forsolvingphysicsproblems.Studentswillbeexpectedtobefamiliarwithbasicprogramming:successful
completionofthe1styearcomputingLab.issufficient. Thereisnorequirementthatthepracticalworkbe
doneusingMicrosoftC++onthedepartmentalcomputers,butanyonewishingtousesomeotherprogram-
minglanguageorcomputershouldconsultthelecturerbeforehand. Thisistomakesurethereisbothhelp
availablefromdemonstratorsandthatitwillbepossibletoassesstheworksatisfactorily.
Adetailedtimetableforthecourseisgivenbelow.
Eachofthe4(5)sectionsofthecoursewillbeaccompaniedbyashortproblemsheetincludingpractical
exercisesandlongerprojects.
Studentswillbeexpectedtodothefirstprojectand2others.Therewillbeachoiceofprojectsavailable
inthelatersections.
Thisyearyouwillberequiredtosubmityourprogramelectronicallyinadditiontoincludingitasan
appendix in your project reports. You should also include an example input file and the corresponding
output file. This will enable the assessors to check and run your program if necessary. Your program
shouldbesentasaplaintexte-mailattachment(s)1toComputational-Physics, A. MacKinnon
[e-mail]
first project is to be handed in to Martin Morris in the Computing Lab by 5.00pm on the 28th October.
(cid:1)(cid:3)(cid:2)(cid:5)(cid:4)(cid:7)(cid:6)
Thesewillbeassessedandreturnedbythedemonstratorswhoshouldprovidethestudentswithfeedback
whichwillbeofuseinthelongerprojectslateron. Thisprojectwillcountfor ofthefinalresult.
1UsingOutlookclick Insert then File
1
During the Lab. sessions staff and demonstrators will be available, but the students are expected to
do the practical work and the problems outside these times as well. I am aware that there is an General
Relativitylecturescheduledfor3.00pmonTuesdays.NeverthelesstheLab.willbereservedforthiscourse
forthewholeafternoonanddemonstratorswillbeavailableduringthisperiod.
Thereportonthe2ndproject,whichshouldbechosenfromthe2offeredonPartialDifferentialEqua-
tions,shouldbesubmittedby5.00pmonNovember18th.
The3rdprojectmaybechosenfromthoseofferedonMatrixandMonteCarlomethods(orpossiblyon
ComputerAlgebra)andthereportshouldbesubmittedby5.00pmonFridayDecember13th,thelastday
ofterm. This timetablehas beendesignedtomakesurethatstudentsareworkingatasteadyrateduring
theterm.
(cid:1)(cid:9)(cid:8)(cid:10)(cid:4)(cid:7)(cid:6)
From10.00am–12.00noononMonday6thJanuary2003,therewillbeashortwrittentestunderexam
(cid:1)(cid:11)(cid:8)(cid:10)(cid:4)(cid:12)(cid:6)
conditions. This will count for of the final result. The 2nd and 3rd projects will also count for
.
Thusthedistributionofmarksforthevariouspartsofthecoursewillbeapproximately
Project1 10
Project2 30
Project3 30
Test 30
Week Lectures Practicals Deadlines
Sep. 30 Thu. 9.00am
Fri. 2.00pm
Fri. 4.00pm
Oct. 7 Fri. 4.00pm Tue. 2.00–5.00pm
Oct. 14 Fri. 4.00pm Tue. 2.00–5.00pm
Oct. 21 Fri. 4.00pm Tue. 2.00–5.00pm
Oct. 28 Fri. 4.00pm Tue. 2.00–5.00pm Mon. 5.00pm1stProject
Nov.4 Tue. 2.00–5.00pm
Nov.11 Thu. 9.00am Tue. 2.00–5.00pm
Nov.18 Thu. 9.00am Tue. 2.00–5.00pm Mon. 5.00pm2ndProject
Nov.25 Thu. 9.00am Tue. 2.00–5.00pm
Dec. 2 Thu. 9.00am Tue. 2.00–5.00pm
Dec. 9 Thu. 9.00am Tue. 2.00–5.00pm Fri. 5.00pm3rdProject
Jan. 6 Test: Mon. 10.00am–12.00noon
For those with surnames beginning with A–K the test will be in lecture theatre 2 and for those with
surnamesbeginningwithL–Zinlecturetheatre3.
Therewillbenoexaminationinthesummer.
AcopyofthetimetableaswellascontactdetailsforthedemonstratorsispostedontheExchangeserver.
ThiscanbeaccessedfromOutlookunder//PublicFolders/AllPublicFolders/Physics/ComputationalPhysics.
FurtherinformationisavailableontheWorldWideWebunder
http://www.cmth.ph.imperial.ac.uk/angus/Lectures/compphys/.
Chapter 1
Ordinary Differential Equations
1.1 Types of Differential Equation
In this chapter we will consider the methods of solution of the sorts of ordinary differential equations
(ODEs)whichoccurverycommonlyinphysics. ByODEswemeanequationsinvolvingderivativeswith
respect to a single variable, usually time. Although we will formulate the discussion in terms of linear
ODEsforwhichweknowtheanalyticalsolution,thisissimplytoenableustomakecomparisonsbetween
thenumericalandanalyticalsolutionsanddoesnotimplyanyrestrictiononthesortsofproblemstowhich
themethodscanbeapplied.Inthepracticalworkyouwillencounterexampleswhichdonotfitneatlyinto
thesecategories.
Theworkin thissectionisalso consideredin chapter16ofPressetal. (1992)orchapterII ofPotter
(1973)
Weconsider3b(cid:13)(cid:15)a(cid:14) sicdifferentialequations:
(cid:14)(cid:23)(cid:22) (cid:4) (cid:24) (cid:14)(cid:25)(cid:22)(cid:26)(cid:14)(cid:28)(cid:27)(cid:30)(cid:29) (cid:31)(cid:7)!#"%$ (cid:16)%&
(cid:13)(cid:17)(cid:16)(cid:19)(cid:18)(cid:21)(cid:20) (cid:20)
(cid:13)(cid:15)(cid:14) TheDecayequation (1.1)
$ (cid:14)(cid:23)(cid:22) (cid:4) (cid:24) (cid:14)(cid:25)(cid:22)(cid:26)(cid:14)(cid:28)(cid:27)(cid:30)(cid:29) (cid:31)(cid:7)!#" (cid:16)%&
(cid:13)(cid:17)(cid:16) (cid:20) (cid:18)’(cid:20)
(cid:13)(cid:15)(cid:14) TheGrowthequation (1.2)
(cid:14),(cid:22) (cid:4) (cid:24) (cid:14)(cid:25)(cid:22)(cid:26)(cid:14)(cid:28)(cid:27)(cid:30)(cid:29) (cid:31)(cid:7)!#".- (cid:16)%&
(cid:13)(cid:17)(cid:16)(cid:19)((cid:21))+* )+*
TheOscillationequation (1.3)
whicharerepresentativeofmostmorecomplexcases.
Higherorderdifferentialequationscanbereducedto1storderbyappropriatechoiceofadditionalvari-
ables. Thesimplestsuchchoiceistodefinenewvariablestorepresentallbutthehighestorderderivative.
Forexample,thedampedharmonicoscillato(cid:13)(cid:17)r0 (cid:14)equatio(cid:13)(cid:15)n(cid:14) ,usuallywrittenas
/ (cid:14)(cid:23)(cid:22) (cid:4)
(cid:13)1(cid:16) 0 (cid:18)32 (cid:13)(cid:17)(cid:16)(cid:19)(cid:18)54
(1.4)
(cid:14) (cid:22)(cid:26)(cid:13)(cid:15)(cid:14)17(cid:28)(cid:13)(cid:17)(cid:16)
6
canberewrittenintermsof andvelocity(cid:13) intheformofapairof1storderODEs
6 2 4 (cid:14)1(cid:22) (cid:4)
(cid:13)(cid:17)(cid:16) (cid:18) / 6 (cid:18) /
(cid:13)(cid:15)(cid:14) (1.5)
$ (cid:22) (cid:4)
(cid:13)1(cid:16) 6
(1.6)
8 8
Similarlyany thorderdifferentialequationcanbereducedto 1storderequations.
(cid:14)
Such systems of ODEs can be written6in a very concise notati8on by defining a vector, y say, whose
elementsaretheunknowns,suchas and i(cid:13)n(1.1). AnyODEin unknownscanthenbewritteninthe
generalform " (cid:16)%&;(cid:22) (cid:4)
(cid:13)(cid:17)(cid:16)9(cid:18) :
y
f y (1.7)
3
OrdinaryDifferentialEquations 4
8
(cid:16)
whereyandfare –componentvectors.
Rememberthatthereisnosignificanceintheuseoftheletter intheaboveequations. Thevariableis
notnecessarilytimebutcouldjustaseasilybespace,asin(1.9),orsomeotherphysicalquantity.
Formallywecanwritethesolutionof(1.7)as
"<(cid:16)%&;(cid:22) "<(cid:16) (cid:27) &=$3>3? " "<(cid:16)CBD& (cid:16)CBE&(cid:15)(cid:13)(cid:17)(cid:16)CB
:
y y ?A@ f y (1.8)
(cid:16)%(cid:27)GFH(cid:16)
"<(cid:16)%&
by integrating both sides over the interval . Although (1.8) is formally correct, in practice it is
usuallyimpossibletoevaluatetheintegralontheright–hand–sideasitpresupposesthesolutiony . We
willhavetoemployanapproximation.
(cid:16) (cid:16)I(cid:22) (cid:4)
All differential equations require boundary conditions. Here we will consider cases in which all the
(cid:16)
boundary conditions are defined at a particular value of (e.g. ). For higher order equations th*e
boundaryconditionsmaybedefinedatdifferentvaluesof . Themodesofaviolinstringatfrequency
(cid:13)(cid:17)0J(cid:14) 0
obeytheequation (cid:22)L$ * (cid:14)
(cid:13)1K 0 M 0
(cid:14)N(cid:22) (cid:4) (1.9)
withboundaryconditionssuchthat atbothendsofthestring. We shallconsidersuchproblemsin
chapter3.
1.2 Euler Method
(cid:16)P(cid:22)(cid:26)(cid:16)RQ(cid:28)S=TU$V(cid:16)CQ
O
Consideranapproximatesolutionof(1.8)overasmallinterval bywritingtheintegralas
> ?AWYX(cid:12)Z " "<(cid:16) B& (cid:16) B&(cid:17)(cid:13)(cid:17)(cid:16) B (cid:1) (cid:16) " "<(cid:16) Q & (cid:16) Q &\[
: O :
? W f y f y (1.10)
"](cid:16)CQ(cid:28)S#TJ&U(cid:22) "](cid:16)CQ1&=$ (cid:16) " "<(cid:16)CQ(cid:17)& (cid:16)CQ1&^[
toobtain O :
y y f y
or,inamoreconcisenotation, Q(cid:10)S#T (cid:22) Q $ (cid:16) " Q (cid:16) Q &;(cid:22) Q $ (cid:16) Q [
O : O
y y f y y f (cid:16) (1.11)
O
Wecanintegrateoveranylargerintervalbysubdividingtherangeintosectionsofwidth andrepeating
(1.11)foreachpart.
Equivalentlywecanconsiderthatweha(cid:13)veapproxQ(cid:28)iS=mTat$edQthederivativewithaforwarddifference
(cid:1) [
(cid:13)(cid:17)(cid:16)=_Q (cid:16)
y_ y O y
_ (1.12)
_
WewillalsocomeacrosscentredandbackwardQ(cid:28)S=diTff$ereQ(cid:7)nce(cid:30)eTs,
(cid:13) c f (cid:16)
(cid:13)(cid:17)(cid:16) _Q (cid:1)a‘b yQ $ OQ(cid:7)egyT centred
y_
_ (cid:16) (1.13)
_ bd y O y
backward
respectively.
Q (cid:22) "<(cid:16)CQ(cid:17)& (cid:16) (cid:16)%Q(cid:23)(cid:22) (cid:16) Q (cid:22) " Q (cid:16)CQ(cid:17)&
Herewehaveusedanotationwhichisveryc8goOmmonincomput:ationalphysics,inwhichwecalculate
y y atdiscretevaluesof givenby ,andf f y .
Inwhatfollowswewilldropthevectornotationyexceptwhenitisimportantforthediscussion.
OrdinaryDifferentialEquations 5
1.2.1 OrderofAccuracy (cid:14)(cid:30)"<(cid:16)%& (cid:16) Q
HowaccurateistheEulermethod? Toquantifyth(cid:13)(cid:17)i(cid:14)swecon(cid:16)s0id(cid:13)er0 (cid:14)aTaylorexpansionof around
(cid:14) Q(cid:28)S=T (cid:22)h(cid:14) Q (cid:16) O [ [ [
(cid:18) O (cid:13)(cid:17)(cid:16) _Q (cid:18) f (cid:13)(cid:17)(cid:16) 0 _Q (cid:18)
_ _
_ _ (1.14)
_ _
andsubstitutethisinto(1.11) (cid:13)(cid:15)(cid:14) (cid:16)C0 (cid:13)(cid:17)0(cid:5)(cid:14)
(cid:14) Q (cid:16) O [ [ [ (cid:1) (cid:14) Q $ (cid:16)%jk"A(cid:14) Q (cid:16) Q &
(cid:18) O (cid:13)(cid:17)(cid:16)i_Q (cid:18) f (cid:13)(cid:17)(cid:16) 0 _Q (cid:18) O :
_ _
_ _ (cid:13)(cid:15)(cid:14) (1.15)
_ _
(cid:22)l(cid:14)(cid:10)Q (cid:16)
(cid:18) O (cid:13)(cid:17)(cid:16) _Q :
_
_ (1.16)
_ (cid:16)
O
wherewehaveused(1.7)toobtainthefinalform. Hence,weseethatthetermin intheexpansionhas
beencorrectlyreproducedbytheapproximation,butthatthehigherorderQ termsarewrong. Wetherefore
(cid:16)
describetheEulermethodas1storder8accurate. O
Anapproximationtoaquantityis thorderaccurateifthetermin intheTaylorexpansionofthe
quantityis correctlyreproduced. The orderofaccuracyofa method is theorder ofaccuracywith which
theunknownisapproximated.
Notethatthetermaccuracyhasaslightlydifferentmeaninginthiscontextfromthatwhichyoumight
use to describethe resultsofan experiment. Sometimesthe term orderof accuracyisused to avoidany
ambiguity.
(cid:16)%0
TheOleading order deviation is called the truncation error. Thus in (1.2.1) the truncation error is the
termin .
1.2.2 Stability
TheEulermethodis1storderaccurate. Howeverthereisanotherimportantconsiderationinanalysingthe
method: stability. Let us suppose that at some time the actual numerical solution deviates from the true
(cid:14)
solutionOofthedifferenceequation(1.11)(N.B.nottheoriginaldifferentialequation(1.7))bysomesmall
amount ,due,forexample,tothefiniteaccuracyofthecomputer.Thejnaddingthisinto(1.11)gives
(cid:14)(cid:10)Q(cid:28)S#T (cid:14)(cid:10)Q(cid:28)S#Tm(cid:22)(cid:26)(cid:14)(cid:10)Q (cid:14)nQo$ (cid:16)qp<jk"<(cid:14)(cid:10)Q (cid:16)CQ(cid:17)& (cid:14)(cid:10)Qnt
(cid:18) O (cid:18) O O : (cid:18)sr (cid:14) _Q O :
_
r _ (1.17)
_
jk"A(cid:14) (cid:16)%& (cid:14)
uv :
(cid:14)
wheretheterminbrackets,O , istheTaylorexpansionof j withrespectto . Subtracting(1.11)we
obtainalinearequationfor (cid:14) Q(cid:28)S#T (cid:22) p (cid:2) $ (cid:16) t (cid:14) Q
O O r (cid:14)w_Q O :
_
r _ (1.18)
_
whichitisconvenienttowriteintheform (cid:14)(cid:10)Q(cid:10)S#Tm(cid:22)lx (cid:14)(cid:10)Qy[
O O
x (cid:14)(cid:7)Q (1.19)
O 8
x (cid:2)
If has a magnitude greater than one then will tend to grow with increaz siz(cid:7)n{g and may eventually
dominateovertherequiredsolution. HencetheEulermejthodisstableonlyif or
$ (cid:2) (cid:2) $ (cid:16) (cid:2) [
{ O r (cid:14) { (cid:18)
r (1.20)
(cid:16)
O
(cid:16)
As ispositivebydefinitionthe2O ndinequalityimpliesthatthederivativemustalsobepositive. The1st
inequalityleadstoarestrictionon ,namely j
(cid:16) f 7 [
O { r (cid:14)
r (1.21)
OrdinaryDifferentialEquations 6
x
z z
x 0 (cid:2)
When thederivativeis complexmoz rez c{are is requiredin the calculation of . Inthiscaseit is easierto
lookforsolutionsofthecondition . F(cid:2)orthe(cid:16) 0osc0illati(cid:2)onequation(1.1c)theconditionbecomes
(cid:18) O * {
(cid:16) (1.22)
O *
which is impossible to fulfil for real and . Comparing these result with our 3 types of differential
equations(1.1)wefindthefollowingstabilityconditions
(cid:16) f 7
DO ec{ ay (cid:20) Growth Oscillation
unstable unstable
TheEulermethodisconditionallystableforthedecayequation.
Amethodisstableifasmalldeviationfrom thetruesolutiondoesnottendtogrowas thesolutionis
iterated.
1.2.3 TheGrowthEquation
Actually,ouranalysisdoesn’tmaketoomuchsenseinthecaseofthegrowthequationasthetruesolution
(cid:14) Q(cid:15)|\Q (cid:14) Q |JQ
shouldgrowanyway. Amore sensibleconditionwouO ld bethatthe relativeerrorin thesolution doesnot
grow. Thiscanbeachievedbysubstituting for aboveandlookingfortheconditionthat does
notgrow.Wewillnottreatthiscasefurtherherebutitis,infact,veryimportantinproblemssuchaschaos,
inwhichsmallchangesintheinitialconditionsleadtosolutionswhichdivergefromoneanother.
1.2.4 ApplicationtoNon–LinearDifferential Equations
Thelinear–differentialequationsinphysicscanoftenbesolvedanalyticallywhereasmostnon–linearones
can only be solvednumerically. It is importanttherefore to be able to apply the ideas developedhereto
suchcases. (cid:13)(cid:17)(cid:14)
Considerthesimpleexample (cid:14) 0 (cid:22) (cid:4) [
(cid:13)(cid:17)(cid:16)(cid:19)(cid:18)(cid:21)(cid:20)
jk"A(cid:14) (cid:16)%&,(cid:22) (cid:14)(cid:15)0 jg7 (cid:14)}(cid:22) f (cid:14) (1.23)
: (cid:20) (cid:20)
r r
In this case and which can be substituted into (1.21) to give the stability
(cid:16) (cid:2) 7(cid:17)" (cid:14)(cid:15)&
condition O { (cid:20) :
(cid:14)
(1.24)
(cid:14) (cid:16)
which depends on , unlike the simpler cases. In writing a program to solOve such an equation it may
therefore be necessary to monitor the value of the solution, , and adjust as necessary to maintain
stability.
1.2.5 ApplicationtoVector Equations
O
7 (cid:127)
Alittlemorecareisrequiredwhenyandfare~#ve(cid:127)(cid:129)(cid:128)ctors. Inthiscase yisanarbitraryinfinitesimalvector
r r (cid:128)
andthederivative f yisamatrixFwithcomponentsj
(cid:22)
(cid:127) (cid:128) r (cid:14)
r (1.25)
j (cid:14) (cid:127)
(cid:127) (cid:127) (cid:128)
(cid:128)
inwhich and representthecompQ(cid:28)oS#nT(cid:132)e(cid:131) ntsoffQ(cid:28)a(cid:131) ndyre(cid:128)spectjively.HeQ(cid:28)n(cid:131) ce(1.18)takestheform
(cid:14)y(cid:130) (cid:22) (cid:14)y(cid:130) $ (cid:16)(cid:17)(cid:133) (cid:14)y(cid:130)
O O O r (cid:14) _Q O
_
r _ (1.26)
Q(cid:28)S#T (cid:22) $ (cid:16) Q (cid:22) _ Q [
O u O v(cid:134)O O
y I F y G y (1.27)
This leads directly to the stability condition that all the eigenvalues of G must have modulus less than
unity(seeproblem6).
Ingeneralanyofthestabilityconditionsderivedinthiscourseforscalarequationscanbere–expressed
inaformsuitableforvectorequationsbyapplyingittoalltheeigenvaluesofanappropriatematrix.
OrdinaryDifferentialEquations 7
1.3 The Leap–Frog Method
How can we improve on the Euler method? The most obvious way would be to replace the forward
differencein(1.12)withacentreddifferQ(cid:10)enS#cT e(cid:22)(1.1Q(cid:7)3e(cid:30))Tto$ gfet(cid:16)th"efQor(cid:16)%m&\[ula
O :
y y f y (1.28)
Q(cid:28)S#T Q(cid:7)egT
Ifweexpandbothy andy (cid:13)(cid:17)(cid:14) asinsec(cid:16)C0tio(cid:13)1n0J1(cid:14) .2.1(1.2(cid:16)C8(cid:135) )(cid:13)(cid:17)b(cid:135)(cid:5)e(cid:14)comes
(cid:14)(cid:10)Q (cid:16) O O [(cid:5)[ [
(cid:18) O (cid:13)(cid:17)(cid:16)i_Q (cid:18) f (cid:13)(cid:17)(cid:16) 0 _Q (cid:18) (cid:136) (cid:13)(cid:17)(cid:16) (cid:135) _Q (cid:18)
_ _ _
_ _ _
(cid:13)(cid:15)_(cid:14) (cid:16)C0 (cid:13)(cid:17)0 _(cid:14) (cid:16)C(cid:135) (cid:13)1(cid:135)J_(cid:14)
(cid:22)(cid:26)(cid:14)(cid:10)Qo$ (cid:16) O $ O [ [ [(cid:12)$ f (cid:16)%j(cid:137)Q
O (cid:13)1(cid:16)P_Q (cid:18) f (cid:13)1(cid:16) 0 _Q (cid:136) (cid:13)(cid:17)(cid:16) (cid:135) _Q (cid:18) O
_ _ _
(cid:13)(cid:15)(cid:14) __ (cid:16)C0 (cid:13)(cid:17)0 (cid:14) __ (cid:16)C(cid:135) (cid:13)1(cid:135)J(cid:14) __ (1.29)
(cid:22)(cid:26)(cid:14)(cid:10)Q (cid:16) O $ O [ [ [
(cid:18) O (cid:13)1(cid:16)P_Q (cid:18) f (cid:13)1(cid:16) 0 _Q (cid:136) (cid:13)(cid:17)(cid:16) (cid:135) _Q (cid:18)
_ _ _
_ _ _ (1.30)
_ _ _
(cid:16) 0
O
(cid:14)(cid:12)Q (cid:14)(cid:7)Q (cid:14)(cid:10)Q(cid:7)egT
from which allterms up to cancelso that the method is clearly2nd order accurate. Note in passing
(cid:14)(cid:12)Q(cid:28)S#T
that using(1.28) 2 consecutivevaluesof are requiredin orderto calculatethenextone: and
arerequiredtocalculate . Hence2boundaryconditionsarerequired,eventhough(1.28)wasderived
from a 1st order differential equation. This so–called leap–frog method is more accurate than the Euler
(cid:14)nQ
method, but is it stable? ROepeating the same analysis (section 1.2.2) as for the Euler method we again
obtainalinearequationfor j
(cid:14)(cid:10)Q(cid:10)S#Tm(cid:22) (cid:14)nQ(cid:7)e(cid:30)TU$ f (cid:16) (cid:14)(cid:10)Qy[
O O O r (cid:14) _Q O
_
r _ (1.31)
_
(cid:14)nQ(cid:23)(cid:22)(cid:138)x (cid:14)(cid:10)Q(cid:7)egT (cid:14)nQ(cid:28)S#Tm(cid:22)hx(cid:7)0 (cid:14)(cid:10)Q(cid:7)e(cid:30)T
O O O O
Weanalysethisequationbywriting and j toobtain
x 0 (cid:22) (cid:2) $ f (cid:16) x
O r (cid:14)w_Q
_
r _ (1.32)
_
j j 0
whichhasthesolutions x,(cid:22) (cid:16) (cid:16) (cid:2) [
O r (cid:14) ((cid:26)(cid:139) (cid:140) O r (cid:14);(cid:141) (cid:18)
r r (1.33)
$ (cid:2)
(cid:2)
The product of the 2 solutions is equal to the constant in th(cid:142)e quadratic equation, i.e. . Since the 2
x (cid:2)
solutions are different, one of them always has magnitude .z Szwin(cid:142) ce for a small random error it is
impossible to guarantee that there will be no contribution with this contribution will tend to
dominateastheequationisiterated.Hencethemethodisunstable.
Thereisanimportantexceptiontothisinstability: whenthepartialderivativeispurelyimaginary(but
x jg7 (cid:14)(cid:143)(cid:22)
notwhenithassomegeneralcomplexvalue),thequantityunderthesquarerootin(1.33)canbe((cid:144)ne)(cid:145)*gative
r r
and both ’s have modulus unity. Hence, for the case of oscillation (1.1c) where , the
(cid:16) (cid:2) 7 [
algorithmisjuststable,aslongas O { *
(1.34)
Thestabilitypropertiesoftheleap–frogmethodaresummarisedbelow
(cid:16) (cid:2) 7
Decay Growth OO sc{ illati*on
unstable unstable
Againthegrowthequationshouldbeanalysedsomewhatdifferently.
OrdinaryDifferentialEquations 8
1.4 The Runge–Kutta Method
So far we have found one method which is stable for the decay equation and another for the oscillatory
equation.Canwecombinetheadvantagesofboth?
T
Asapossiblecompromisecons(cid:14)(cid:15)iQ(cid:10)BdS#erT(cid:132)(cid:146)the(cid:22)(cid:26)fo(cid:14)(cid:10)lQolow$ ing(cid:16)%tjkw"Ao(cid:14)(cid:10)Qste(cid:16)CpQ(cid:17)a& lgorithm(ignoringvectors)
0 0 O :
(cid:14)(cid:10)Q(cid:28)S#Tm(cid:22)(cid:26)(cid:14)(cid:10)Qo$ O (cid:16)%jk"<(cid:14)(cid:7)Q(cid:28)B S#TR(cid:146) 0 : (cid:16) Q(cid:10)S#T(cid:132)(cid:146) 0 &\[ (1.35)
(cid:14) Q(cid:28)B S#TR(cid:146) (1.36)
0
Inpracticetheintermediatevalue isdiscardedaftereachstep. Weseethatthismethodconsistsof
anEuler(seesection1.2)stepfollowedbyaLeap–Frog(seesection1.3)step. Thisiscalledthe2ndorder
Runge–Kuttaortwo–stepmethod.Itisinfactoneofahierarchyofrelatedmethodsofdifferentaccuracies.
Thestabilityanalysisfor(1.4)iscarriedoutinthj esamT ewayjasb0\e(cid:149) fore.Herewesimplyquotetheresult
(cid:14) Q(cid:28)S#T (cid:22)(cid:148)(cid:147) (cid:2) $ (cid:16) (cid:16) (cid:14) Q [
O O r (cid:14) (cid:18) 0 (cid:140) O r (cid:14);(cid:141) O
r r (1.37)
jg7 (cid:14) (cid:16)
r r
Inderivingthisresultitisnecessarytoassumethatthederivatives, areindependentof . Thisisnot
usuallyaproblem.
From(1.37)weconcludethestabilityconditions
T
(cid:16) f 7 (cid:2) " (cid:16) & (cid:150) (cid:2)
DO ec{ ay (cid:20) Growth Os(cid:18) cil(cid:150)laOtio*n {
unstable
(cid:16)(cid:23)(cid:151) (cid:2)
Notethatintheoscillatorycasethemeth*oOdisstrictlyspeakingunstablebuttheeffectissosmallthat
it can be ignored in most cases, as long as . This method is often used for damped oscillatory
equations.
1.5 The Predictor–Corrector Method
This method is very similar to and often confused with the Runge–Kutta (see section 1.4) method. We
T
considersubstitutingthetrapezQ(cid:10)oS#idTal(cid:22) rulQef$orth(cid:16);e(cid:152)es"tiQ(cid:28)mS#aTte(cid:16)CoQ(cid:28)fS#tT h&einte" gQral(cid:16)CQ1in&C(cid:153)#(1[ .8)toobtaintheequation
0(cid:12)O : (cid:18) :
y y f y f y (1.38)
Q(cid:28)S#T
Unfortunately the presence of y on the right hand side makes a direct solution of (1.38) impossible
exceptforspecialcases(sees(cid:14)eQ(cid:10)BctS#ioT n(cid:22)h1.(cid:14)n6Q(cid:19))b$ elo(cid:16)%wjk."A(cid:14)(cid:10)AQ p(cid:16)CoQ(cid:17)s&siblecompromiseisthefollowingmethod
OT :
(cid:14) Q(cid:10)S#T (cid:22)h(cid:14) Q $ (cid:16) (cid:152)jk"<(cid:14) Q(cid:28)B S#T (cid:16) Q(cid:28)S#T jk"A(cid:14) Q (cid:16) Q & (cid:153) [
0 O : (cid:18) : (1.39)
(cid:14)(cid:12)Q(cid:28)S=T (1.40)
This method consists of a guess for based on the Euler (see section 1.2) method (the Prediction)
followedbyacorrectionusingthetrapezoidalruletosolvetheintegralequation(1.8). Theaccuracyand
stabilitypropertiesareidenticaltothoseoftheRunge–Kuttamethod.
1.6 The Intrinsic Method
Returning to the possibility of solving the integral equation(1.8) using the trapezoidal or trapezium rule
(1.38),letusconsiderthecaseofalineardifferentiTalequation,suchasourexamples(1.1). Forthedecay
(cid:14)(cid:10)Q(cid:28)S=Tm(cid:22)h(cid:14)(cid:10)Qo$ (cid:16) (cid:14)(cid:10)Q(cid:10)S#T (cid:14)(cid:10)Q
equationwehave 0 O u(cid:20) (cid:18)(cid:21)(cid:20) v
(1.41)
OrdinaryDifferentialEquations 9
(cid:14)(cid:12)Q(cid:10)S#T (cid:14)nQ
T
whichcanberearrangedintoanexplicitequationf(cid:2) o$r (cid:16)g(cid:154) asafunctionof
(cid:14) Q(cid:28)S#T (cid:22) p 0T O (cid:20) t (cid:14) Q [
(cid:2) (cid:16)g(cid:154)
(cid:18) 0(cid:12)O (cid:20)
(1.42)
x
This intrinsicmethod is2nd order accurate asthat is the accuracyofthe trapezoidalrule for integration.
(cid:151) (cid:2)
Whataboutthestability? Applyinguv thesamemethodologyasbeforewefindthatthecrucialquantity, ,is
F
theexpressioninsquarebrackets, ,in(1.42)w(cid:20)hhic(cid:155) hi((cid:144)sa)+*lways forthedecayequationandhasmodulus
unityintheoscillatorycase(aftersubstituting ). Henceitisstableinbothcases. Whyisitnot
usedinsteadoftheothermethods? Unfortunatelyonlyasmallgroupofequations,suchasourexamples,
canberearrangedinthisway.Fornon–linearequationsitmaybeimpossible,andevenforlinearequations
when y is a vector, there may be a formal solution which is not useful in practice. It is always possible
to solve the resulting non–linear equation iteratively, using (e.g. ) Newton–Raphson iteration, but this is
usuallynotworthwhileinpractice.
Infacttheintrinsicmethodisalsostableforthegrowthequationwhenitisanalysedasdiscussedearlier
(section1.2.3),sothatthemethodisinfactstableforall3classesofequations.
Decay Growth Oscillation
stable stable stable
1.7 Summary
Inthesenoteswehaveintroducedsomemethodsforsolvingordinarydifferentialequations. However,by
farthemostimportantlessontobelearnedisthattobeusefulamethodmustbebothaccurateandstable.
The latter condition is often the most difficult to fulfil and careful analye(cid:30)sTis of the equations to be solved
(cid:16) (cid:16)
maybenecessarybeforechoosinganappropriatemethodtousefOor{fin*dinganumericalsolution. O
Thestabilityconditionsderivedabovetendtohavetheform whichmaybeinterpretedas
(cid:16)
shouldbelessthanthecharacteristicperiodofoscillation.OThisconformswithcommonsense. Infactwe
canwritedownamoregeneralcommonsensecondition: shouldbesmallcomparedwiththesmallest
timescalepresentintheproblem.
Finally, many realistic problems don’tfall into the neat categories of (1.1). The simplest example is
a damped harmonic oscillator. Often it is difficult to find an exact analytical solution for the stability
condition. Itpaysinsuchcasestoconsidersomeextremeconditions,suchas(e.g.) veryweakdamping
(cid:14)
or very strong damping, work out the conditions for these cases and simply choose the most stringent
condition. Innon–linearproblemsthecaseswhentheunknown, ,isverylargeorverysmallmayprovide
tractablesolutions.
1.8 Problems
1. Write down definitions of the terms order of accuracy, truncation error, conditional stability as
appliedtothenumericalsolutionofordinarydifferentialequations.
2. Writedown1storderaccurate(cid:13)1fijnitediffer(cid:13)(cid:17)e0Ynjceapprox(cid:13)(cid:17)im(cid:135)(cid:5)jationsfor j
(cid:13)(cid:17)K (cid:13)(cid:17)K 0 (cid:13)(cid:17)K (cid:135) r (cid:16) r K
a) b) c) d) r r
(cid:14) Q(cid:28)S=T (cid:14) Q M (cid:14) Q(cid:15)e(cid:30)T (cid:14) (cid:14) Q
(cid:156) (cid:18)(cid:158)(cid:157) (cid:18)
Hint:theresulthastobesomethinglike .Expandthe saround andchoose
thecoefficientstoeliminatecontributionsfromunwantedtermsintheexpansion.
j
N.B.Thisquestionreferstotheaccuracyoftheapproximationforthederivativesgiven,nottothe
accuracyof .