MacKinnon - Computational Physics (523159), страница 2
Текст из файла (страница 2)
Hence, we see that the term in O in the expansion hasQbeen correctly reproduced by the approximation, but that the higher order terms are wrong. We thereforedescribe the Euler method as 1st order accurate.An approximation to a quantity is th order accurate if the term inin the Taylor expansion of thequantity is correctly reproduced. The order of accuracy of a method is the order of accuracy with whichthe unknown is approximated.Note that the term accuracy has a slightly different meaning in this context from that which you mightuse to describe the results of an experiment.
Sometimes the term order of accuracy is used to avoid anyambiguity.The leading order deviation is called the truncation error. Thus in (1.2.1) the truncation error is theterm in.8OO %01.2.2 StabilityThe Euler method is 1st order accurate. However there is another important consideration in analysing themethod: stability.
Let us suppose that at some time the actual numerical solution deviates from the truesolution of the difference equation (1.11) (N.B. not the original differential equation (1.7)) by some smallamount , due, for example, to the finite accuracy of the computer. Then adding this into (1.11) givesOQS#T O QS#TmQ O nQo$ O qp<jk"<Q : CQ& sr j _ Q O Qnt :(1.17)__r _ kjA"%&: with respect to . Subtracting (1.11) wewhere the term in brackets, u v , is the Taylor expansion ofobtain a linear equation for OjO QS#T p $ O r w__ Q t O Q :(1.18)__rwhich it is convenient to write in the form(1.19)O QS#Tmlx O Qy[xQ will tend to grow with increasing 8 and may eventuallyIf has a magnitude greater than one then Oxdominate over the required solution. Hence the Euler method is stable only if z z{ or$ { $ O r j { [(1.20)rAs O is positive by definition the 2nd inequality implies that the derivative must also be positive.
The 1stinequality leads to a restriction on O , namelyO { f 7 r j [(1.21)rOrdinary Differential Equations6zx zxz z 0 { O 0* 0 { (1.22)which is impossible to fulfil for real O and * . Comparing these result with our 3 types of differentialWhen the derivative is complex more care is required in the calculation of .
In this case it is easier tolook for solutions of the condition. For the oscillation equation (1.1c) the condition becomesequations (1.1) we find the following stability conditionsO { f7DecayGrowthunstableOscillationunstableThe Euler method is conditionally stable for the decay equation.A method is stable if a small deviation from the true solution does not tend to grow as the solution isiterated.1.2.3 The Growth EquationActually, our analysis doesn’t make too much sense in the case of the growth equation as the true solutionshould grow anyway.
A more sensible condition would be that the relative error in the solution does notgrow. This can be achieved by substitutingforabove and looking for the condition thatdoesnot grow. We will not treat this case further here but it is, in fact, very important in problems such as chaos,in which small changes in the initial conditions lead to solutions which diverge from one another. Q|\QOQ|JQ1.2.4 Application to Non–Linear Differential EquationsThe linear–differential equations in physics can often be solved analytically whereas most non–linear onescan only be solved numerically.
It is important therefore to be able to apply the ideas developed here tosuch cases.Consider the simple example 0 [(1.23)kjA"%,&0gj7}f: and r r which can be substituted into (1.21) to give the stabilityIn this casecondition { 7" & :O(1.24)which depends on , unlike the simpler cases. In writing a program to solve such an equation it maytherefore be necessary to monitor the value of the solution, , and adjust O as necessary to maintainstability.1.2.5 Application to Vector Equations7~#OA little more care is required when y and f are vectors. In this case y is an arbitrary infinitesimal vectorand the derivative f y is a matrix F with componentsr rjr rj (1.18) takes the formin which and represent the componentsy respectively.
Hence QS#T of fQandO y O y $ O r j _ Q O y Q__r$O yQS#T u I O FvO yQ _ GO yQ [(1.25)(1.26)(1.27)This leads directly to the stability condition that all the eigenvalues of G must have modulus less thanunity (see problem 6).In general any of the stability conditions derived in this course for scalar equations can be re–expressedin a form suitable for vector equations by applying it to all the eigenvalues of an appropriate matrix.Ordinary Differential Equations71.3 The Leap–Frog MethodHow can we improve on the Euler method? The most obvious way would be to replace the forwarddifference in (1.12) with a centred difference (1.13) to get the formulayQS#T yQeT $ f O f " yQ : %&\[(1.28)QS#T and yQegT as in section 1.2.1 (1.28) becomesQ O i_ Q O f C0 10J 0 _ Q O C _ Q [[ [__ __ __ C0C01J__Qo$ O 1P_ Q O f 1 0 _ Q $ O _ _ Q [ [ [$ f O %jQ(1.29)______Q O 1 P__ Q O f C0 10 0 __ Q $ O C 1J __ Q [ [ [(1.30)________ 0 cancel_ so that the methodfrom which all terms up to Oisclearly2nd order accurate.
NoteQ are required in order to calculateQ inandpassingQegTthat using (1.28) 2 consecutive values ofthe next one:Q#STare required to calculate. Hence 2 boundary conditions are required, even though (1.28) was derivedIf we expand both yfrom a 1st order differential equation.
This so–called leap–frog method is more accurate than the Eulermethod, but is it stable? Repeating the same analysis (section 1.2.2) as for the Euler method we againobtain a linear equation forO nQO QS#Tm O nQeTU$ f O r j _ Q O Qy[(1.31)__rnQx O QegT and O nQS#Tmh_ x0 O QeT to obtainWe analyse this equation by writing Ox 0 $ f O r jw_ Q x(1.32)r ___which has the solutionsx, O r j ( O r j; 0 [(1.33)rr$The product of the 2 solutions is equal to the constant in the quadratic equation, i.e.
. Since the 2solutions are different, one of them always has magnitude . x Since for a small random error it isimpossible to guarantee that there will be no contribution with z zw this contribution will tend todominate as the equation is iterated. Hence the method is unstable.There is an important exception to this instability: when the partial derivative is purely imaginary (butnot when it has some general complex value), the quantity under the square root in (1.33) can be negativeand both ’s have modulus unity. Hence, for the case of oscillation (1.1c) where, thealgorithm is just stable, as long as(1.34)xjg7 ()*r rO { 7* [The stability properties of the leap–frog method are summarised belowDecayunstableGrowthunstableO { 7*OscillationAgain the growth equation should be analysed somewhat differently.Ordinary Differential Equations81.4 The Runge–Kutta MethodSo far we have found one method which is stable for the decay equation and another for the oscillatoryequation.
Can we combine the advantages of both?As a possible compromise consider the following two step algorithm (ignoring vectors)QB S#T 0 Qo$ 0T O %jk"AQ : CQ&(1.35)QS#TmQo$ O %jk"<QB S#TR 0 : QS#T 0 &\[(1.36) QB S#TR 0 is discarded after each step.
We see that this method consists ofIn practice the intermediate valuean Euler (see section 1.2) step followed by a Leap–Frog (see section 1.3) step. This is called the 2nd orderRunge–Kutta or two–step method. It is in fact one of a hierarchy of related methods of different accuracies.The stability analysis for (1.4) is carried out in the same way as before. Here we simply quote the resultj T j \0 O QS#T $ O r 0 O r ; O Q [(1.37)rr jg7 In deriving this result it is necessary to assume that the derivatives,are independent of .
This is notr rusually a problem.From (1.37) we conclude the stability conditionsO { f7DecayGrowthunstable* O T " O * & { OscillationNote that in the oscillatory case the method is strictly speaking unstable but the effect is so small thatit can be ignored in most cases, as long as. This method is often used for damped oscillatoryequations.1.5 The Predictor–Corrector MethodThis method is very similar to and often confused with the Runge–Kutta (see section 1.4) method. Weconsider substituting the trapezoidal rule for the estimate of the integral in (1.8) to obtain the equationQS#T yQ $ 0T O ; f " yQS#T : CQS#T & f " yQ : CQ1&C#[(1.38)QS#T on the right hand side makes a direct solution of (1.38) impossibleUnfortunately the presence of yyexcept for special cases (see section 1.6) below. A possible compromise is the following method QB S#T h nQ$ O T %jk"AQ : CQ&(1.39) QS#T h Q $ 0 O jk"< QB S#T : QS#T jk"A Q : Q & [(1.40)Q=STThis method consists of a guess forbased on the Euler (see section 1.2) method (the Prediction)followed by a correction using the trapezoidal rule to solve the integral equation (1.8).
The accuracy andstability properties are identical to those of the Runge–Kutta method.1.6 The Intrinsic MethodReturning to the possibility of solving the integral equation (1.8) using the trapezoidal or trapezium rule(1.38), let us consider the case of a linear differential equation, such as our examples (1.1).
For the decayequation we have(1.41)QS=TmhQo$ 0T O u QS#T Q vOrdinary Differential Equations9which can be rearranged into an explicit equation for QS#T p $QS#T as a function of nQT0T O gg t Q [0O (1.42)xThis intrinsic method is 2nd order accurate as that is the accuracy of the trapezoidal rule for integration.What about the stability? Applying the same methodology as before we find that the crucial quantity, , isthe expression in square brackets, , in (1.42) which is alwaysfor the decay equation and has modulusunity in the oscillatory case (after substituting).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
