MacKinnon - Computational Physics (523159), страница 4
Текст из файла (страница 4)
The secondequation (2.3b) can then be solved by using the known solution of the first. Note that we could equallyhave chosen the equations the other way round, with the signs of the velocity interchanged.As the 2 equations (2.3) are so similar we expect the stability properties to be the same. We thereforeconcentrate on (2.3a) which is known as the Advective equation and is in fact the conservation of massequation of an incompressible fluid (section 2.6.1)ryÆ M ryÆK [r rNote also that the boundary conditions will usually be specified in the form"<K : %&; "<Ky& at k [ÆÆ<"K%&Kwhich gives the value ofÆ : for all at a particular time.M(2.9)(2.10)2.3.1 A Simple AlgorithmQ QM $ f O K3Ç S#T $ eT^È(2.11)ÆÆ O Æ Æ Qwhere the subscripts represent space steps and the superscripts time steps.
By analogy with the discussion of the Euler (section 1.2) and Leap–Frog K (section 1.3) methods we can see immediately that thismethod is 1st order accurate in and 2nd order in . We note firstly thatQS#T QQQ$ O ryÆ _ 0T O 0 r 0 Æ0 _ [ [ [(2.12)____ÆÆrr_Q_ QQ QwhereasS#T $ egT f O K ryÆK _ T O K r K Æ _ [ [ [(2.13)__Æ Æ__r _r _QS#TQAs a first attempt to solve (2.9) we consider using centred differences for the space derivative and Euler’smethod for the time partPartial Differential Equations15QQQQO yr Æ _ 0T O 0 r 0 Æ0 _ [ [ [ f M O KÊ É f O K ryÆK _ T O K r K Æ _ [[ [(2.14)Or ___r ___r ___r ___so that, when the original differential equation (2.9) is subtracted, we are left with a truncation error whichis 2nd order in the time but 3rd order in the spatial part.
, integration rather than theQ space,Q K . We analyse this byThe stability is a property of the time,K 6 !#" )DË K & to obtainconsidering a plane wave solution for the –dependence by substitutingÆQS=TQM Q6 ¿YÌÎÍ ª Ï 6 ¿YÌÎÍ ªÏ $ f O O K 6 Ç ¿ Í ªÏ XZ $ ¿YÌ Í ªÏCÐ Z È(2.15)or, after dividing out the common exponential factor,QS#TQM 6 p $ ) O O K µ )E¹ " Ë O KÑ& t 6 [(2.16)Q lawQ (sectionQ 2.6) the solutionSince the wave (2.5) and advection (2.9) equations express a conservationQ as a function of time. If we substitute 6 F 6 O6 and subtract (2.16)should neither grow nor decaywe obtain an equation for O6QS=T p $ M O µ " Ky& t Q) O K )E¹ Ë O O6(2.17)O6and substitute these forms into (2.11) to obtainwhich is identical to (2.16).
Hence the stability condition is simply given by the quantity in square brackets.We call this the amplification factor and write it asx $ )EwhereM O O K µ )E¹ " Ë O KÑ&\[xAs is complex the stability condition becomesz x z 0 0 { for all Ë [The condition must be fulfilled for all wave vectors Ë ; otherwise a component with a particular Ë(2.18)(2.19)(2.20)will tendto grow at the expense of the others.
This is known as the von Neumann stability condition 2. Unfortunatelyit is never fulfilled for the simple method applied to the advection equation: i.e. the method is unstable forthe advection equation.2.3.2 An Improved Algorithm — the Lax methodOKOWe see from (2.20) that our simple method is in fact unstable for the advection equation, for all finite valuesofand .
How might we improve on this? Let us consider a minor (?) modification of (2.11)ÆQ QQS=T T Q QMQ Æ 0 Ç Æ #S T Æ egT^È $ f OO K Ç Æ S#T $ Æ egTJÈ [(2.21)in which the term inhas been replaced by an average over its 2 neighbours. When we apply the same(von Neumann) analysis to this algorithm we findQS#T ]p ³J´µ " KÑ&k$VÒ M O µ " Ky& t QË O O K )E¹ Ë O 66(2.22)2 von Neumann is perhaps better known for his work on the theory of computing, but his name is also associated with various otheraspects of mathematical physicsPartial Differential Equations16xz z 0 J³ ´µ 0 " Ë O KÑ& ¼ M O K O $ µ )E¹ 0 " Ë O KÑ&kÓ $so thatwhich is stable for allËO K NÕ MOas long as0 µ 0 " KÑ&)E¹ Ë OM O 01Ô O Kw(2.23)(2.24)(2.25)which is an example of the Courant–Friedrichs–Lewy condition applicable to hyperbolic equations.There is a simple physical explanation of this condition: if we start with the initial condition such thateverywhere except at one point on the spatial grid, then a pointsteps away on thegrid will remain zero until at least time steps later.
If, however, the equation is supposed to describe aphysical phenomenon which travels faster than that then something must go wrong. This is equivalent tothe condition that the time step, , must be smaller than the time taken for the wave to travel the distanceof the spatial step, ; or that the speed of propagation of information on the grid,, must be greaterthan any other speed in the problem.Æ"<K : , & OOK//O K£7 O FÖ O shouldKWhen applying the von NeumannanalysistoK non–linearequations the substitutionJ ²!"&O )DË . should thenbe performedbefore setting OÒ be treated as a constant, independent of(or ).
The resulting stability condition will then depend on , just as in the case of ODE’s (section 1.2.4).2.3.3 Non–Linear Equations2.3.4 Other methods for Hyperbolic EquationsThere is a large number of algorithms applicable to hyperbolic equations in general. As before, theirsuitability for solving a particular problem must be based on an analysis of their accuracy and stability forany wavelength. Here we simply name a few of them:The Lelevier method.The one–step and two–step Lax–Wendroff methods.The Leap–Frog method.The Quasi–Second Order method.These are discussed in more detail in chapter II of Potter (1973) and chapter 17 of Press et al.
(1992).2.4 Eulerian and Lagrangian MethodsIn the methods discussed so far the differential equations have been discretised by defining the value ofthe unknown at fixed points in space. Such methods are known as Eulerian methods. We have confinedourselves to a regular grid of points, but sometimes it is advantageous to choose some other grid. Anobvious example is when there is some symmetry in the problem, such as cylindrical or spherical: oftenit is appropriate to base our choice of grid on cylindrical or spherical coordinates rather than the Cartesianones used so far.Suppose, however, we are dealing with a problem in electromagnetism or optics, in which the dielectricconstant varies in space. Then it might be appropriate to choose a grid in which the points are more closelyspaced in the regions of higher dielectric constant.
In that way we could take account of the fact that thewavelengths expected in the solution will be shorter. Such an approach is known as an adaptive grid.In fact it is not necessary for the grid to be fixed in space. In fluid mechanics, for example, it is oftenbetter to define a volume of space containing a fixed mass of fluid and to let the boundaries of these cellsPartial Differential Equations17move in response to the dynamics of the fluid. The differential equation is transformed into a form in whichthe variables are the positions of the boundaries of the cells rather than the quantity of fluid in each cell.
Asimple example of such a Lagrangian method is described in the Lagrangian Fluid project (section 2.9).2.5 Parabolic Equations — DiffusionParabolic equations such as the diffusion and time–dependent Schr ödinger equations are similar to thehyperbolic case in that the boundary is open and we consider integration with respect to time, but, as weshall see, they present some extra problems.2.5.1 A Simple Method0 $yr Æ U4 r K Æ 0(2.26)rrand discretise it using the Euler methodQ andQ theQ simplest centred 2nd order derivativeQS#T forQ the time derivativeto obtain 4 K O 0 Ç S=T $ f eT È [(2.27)ÆÆ O ÆÆ ÆKApplying the von Neumann analysis to this system by considering a single Fourier mode in space, weobtainKQS#T Q(2.28)6 6 p $×O 4 K O 0 µ )E¹ 0 Ë Of tso that the condition that the method is stable for all Ë givesO { 0T O K10 [(2.29)4Although the method is in fact conditionally stable the condition (2.29) hides an uncomfortableK property:namely, that if we want to improve accuracy and allow for smaller wavelengths by halving O we mustdivide O by × .
Hence, the number of space steps is doubled and the number of time steps is quadrupled:the time required is multiplied by Ø . Note that this is different from the sorts of conditions we haveWe consider the diffusion equation and apply the same simple method we tried for the hyperbolic case.encountered up to now, in that it doesn’t depend on any real physical time scale of the problem.2.5.2 The Dufort–Frankel MethodWe consider here one of many alternative algorithms which have been designed to overcome the stabilityproblems of the simple algorithm.
The Dufort–Frankel method is a trick which exploits the unconditionalstability of the intrinsic method (section 1.6) for simple differential equations.We modify (2.27) to readQeT f QQS#T QeT Q 4 K O 0lÙ S=T $ Ç È eTÚ [(2.30)ÆÆ QS#T O ÆÆÆÆwhich can be solved explicitly forpointQ QQS#T Æ $ at eachQmesheT S#T egTJÈ(2.31)Æ5 Æ Û Ç Æ Æwhere f 4 K O 0 [(2.32)OWhen the usual von Neumann analysis is appliedK to this method it is found to be unconditionally stable.Note however that this does not imply that O and O can be made indefinitely large; common sense tellsQS#TPartial Differential Equations18us that they must be small compared to any real physical time or length scales in the problem. We muststill worry about the accuracy of the method.
Another difficulty this method shares with the Leap–Frog(section 1.3) method is that it requires boundary conditions at 2 times rather than one, even though theoriginal diffusion equation is only 1st order in time.2.5.3 Other MethodsAnother important method for solving parabolic equations is the Crank–Nicholson method, which we shallnot discuss here but which is considered in chapter II of Potter (1973).2.6 Conservative MethodsMost of the differential equations we encounter in physics embody some sort of conservation law.
Itis important therefore to try to ensure that the method chosen to solve the equation obeys the underlyingconservation law exactly and not just approximately. This principle can also have the side effect of ensuringstability. The simplest approach to deriving a method with such properties is to go back to the originalderivation of the differential equation.2.6.1 The Equation of ContinuityThe archetypal example of a differential equation which implies a conservation law is the equation ofcontinuity, which, in its differential form says thatÜryÜ Ý j (2.33)rrepresents a density and j a current density.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
