MacKinnon - Computational Physics (523159), страница 3
Текст из файла (страница 3)
Hence it is stable in both cases. Why is it notused instead of the other methods? Unfortunately only a small group of equations, such as our examples,can be rearranged in this way. For non–linear equations it may be impossible, and even for linear equationswhen y is a vector, there may be a formal solution which is not useful in practice. It is always possibleto solve the resulting non–linear equation iteratively, using (e.g. ) Newton–Raphson iteration, but this isusually not worthwhile in practice.In fact the intrinsic method is also stable for the growth equation when it is analysed as discussed earlier(section 1.2.3), so that the method is in fact stable for all 3 classes of equations.uv hF ()+*DecaystableGrowthstableOscillationstable1.7 SummaryIn these notes we have introduced some methods for solving ordinary differential equations.
However, byfar the most important lesson to be learned is that to be useful a method must be both accurate and stable.The latter condition is often the most difficult to fulfil and careful analysis of the equations to be solvedmay be necessary before choosing an appropriate method to use for finding a numerical solution.The stability conditions derived above tend to have the formwhich may be interpreted asshould be less than the characteristic period of oscillation. This conforms with common sense. In fact wecan write down a more general common sense condition:should be small compared with the smallesttime scale present in the problem.Finally, many realistic problems don’t fall into the neat categories of (1.1). The simplest example isa damped harmonic oscillator. Often it is difficult to find an exact analytical solution for the stabilitycondition. It pays in such cases to consider some extreme conditions, such as (e.g.
) very weak dampingor very strong damping, work out the conditions for these cases and simply choose the most stringentcondition. In non–linear problems the cases when the unknown, , is very large or very small may providetractable solutions.OO { *eTO1.8 Problems1. Write down definitions of the terms order of accuracy, truncation error, conditional stability asapplied to the numerical solution of ordinary differential equations.2. Write down 1st order accurate finite difference approximations for0YjjjK 0 c) K d) r r Kr r QQ=STQQeTM Hint: the result has to be something like .
Expand the s arounda)1jKb)and choosethe coefficients to eliminate contributions from unwanted terms in the expansion.jN.B. This question refers to the accuracy of the approximation for the derivatives given, not to theaccuracy of .Ordinary Differential Equations1j a) K _ Q0j ___ c) K 0 _ Q__1j _ e) K _ Q___10jQS#T$}jQ1j j QS#TU$}jQegTKb) K _ Qf KjQS=O TU$ f jQ jQegTj ___ jQS 0 O $ f jQS#T f jQeTU$}jQe 0d) K _ QfOK__OK0$ T T 0 jQS 0 0 jQS#T$ 0 jQeT T T 0 jQ_ e 0OK3.
Derive expressions for the truncation error of the following difference approximations. [K }¡4. The torsion of a bar is described by the differential equationShow how to re–express this as a system of first order differential equations. L¢5. Write down an expression for solving the differential equationby Euler’s method and show under what conditions the method is stable.£"]k &Write and test a short program (it should only require a few lines) to test the method. Varyto check the validity of the stability condition you have derived.O ,¢and6. Using equation 1.2.5 show that the Euler method is stable for a vector equation provided all theeigenvalues of G have modulus less than or equal to unity.7.
Show that equation 1.37 gives the correct stability condition for both the Runge–Kutta and Predictor–Corrector methods. Why do you think this is a good method for damped oscillatory equations? (Thelast part doesn’t have to be mathematically rigorous).1.9 Project — Classical Electrons in a Magnetic FieldThe simplest source of ODEs in physics is classical mechanics, so we choose such a problem. The dynamics of a charge in a magnetic field is described by the equation/¤/ v ¤ v ¥ B $ v¢(1.43)where and are the mass and charge of the particle respectively, B is the magnetic field andsome sort of friction.¢represents1.9.1 A Uniform FieldBefore considering the more general problem we start with a particle in a spatially uniform field. Thisproblem is analytically solvable and can be used to test the various methods before applying them to thegeneral case.Units¤¦,7 /¢7/Note firstly that there are only 2 independent constants in the problem,and, and that theseconstants have the units of inverse time; in fact the former is the cyclotron frequency and the latter is adamping rate.
In general in any programming problem it pays to think carefully about the units to be usedin the program. There are several reasons for this.Ordinary Differential Equations11§If inappropriate units are used the program may not work at all. An example of this would be theuse of SI units to study the dynamics of galaxies or to study atomic physics. In the formermighteasily arise and be bigger than the largest number representable on the machine, whereas in the lattermay be smaller than the smallest number on the machine and be set automatically to zero withdisastrous consequences.©¨ The problem often has its own natural units and it makes sense to work in these units. This has theconsequence that most of the numbers in your program will be of order unity rather than very largeor very small.In general you should look for the natural units of a problem and write your program appropriately.
Notethat these will generally not be SI or cgs.and. If we decideIn the problem we are considering here there are 2 natural time scales,to work in one of these, e.g. the cyclotron period, we can rewrite (1.43) in the simpler form 6ª B 6« B/ 7 ¤¦ 6« $$ 6 ª $/ 7¤¦__ ¢¤ ¦____ ¢¤ ¦____ 6 ª____ 6«__/ 7¢(1.44)(1.45) 6ª $ $ B 6 « __ ¢¤ ¦__ 6ª(1.46)__ 6« $ _ _ B 6 ª __ ¢¤ ¦ __ 6«(1.47)__¤¦,7 / . Here ¬ B z / 7¤¦ z _and we_ have chosen our coordinate system suchdepending on the sign ofthat the magnetic field, B, is in the –direction. Note, in addition, that choosing the units appropriately hasor perhapseliminated all but one of the constants from the problem. This cuts down on superfluous arithmetic in theprogram.The Analytical SolutionIn order to understand the behaviour of the various methods for ODEs we need to know the analyticalsolution of the problem.
2 dimensional problems such as this one are often most easily solved by turningthe 2D vector into a complex number. Thus by definingwe can rewrite (1.43) in the form6 ® n6 ª }) 6 « 6 ® ¯ $ ¤ ¦ $1 ) 6 ® / 6® ° /l¢ ±(1.48)which can be easily solved using the integrating factor method to give¤ ¦6 ® 6 ® ² !Np$ ) / #$ ° /l¢ ± ]t'[Finally we take real and imaginary parts to find the 6nª and 6 « components¤¦6ª 6 ²³J´µp / 5¶ Jt !N·]$ ° /l¢ ± .¸¤¦6« L$ 6 µ )E¹ p / 5¶ ^tmJ ! · $ ° /l¢ ± ¸(1.49)(1.50)(1.51)Ordinary Differential Equations12Choosing an AlgorithmThe problem as posed does not fit neatly into the categories defined in (1.1).
By considering the accuracyand stability properties of the various algorithms described in these notes you have to decide which is themost appropriate algorithm to use for solving the problem. The computer time required by the algorithm aswell as the ease with which it can be programmed may be legitimate considerations. It may not be possibleto solve the the stability equations exactly in the most general case. Nevertheless you should be able todeduce enough information on which to base a decision (See (e.g.) section 1.7).In your report you should discuss the merits and demerits of the various possible algorithms and explainthe rationale behind your choice. You should also write a program to implement your chosen algorithmand test it for various values of and in all the physically significant regimes.O1.9.2 Crossed Electric and Magnetic FieldsK"]º ª: : &You are now in a position to apply your chosen algorithm to a more complicated problem.
In addition tothe uniform magnetic field, B, we now add an electric field in the –direction, E. Thus (1.43a)must be modified to read 6 ª ¤¦ / 1 » 6« $ ¤¦ / ¤º / 6 « $ ° / ¢ ± 6 ª ¼6 ª $ ° /¢ ± 6 «(1.52)(1.53)You should now write a program to solve (1.9.2) using the most appropriate method as found in section 1.9.1. Try to investigate the behaviour of the system in various physical regimes.
You should also varyto check whether the stability conforms to your expectations. Think about the physical system you aredescribing and whether your results are consistent with the behaviour you would expect.O1.9.3 Oscillating Electric Fieldº ª ºm ³ ´nµ * Finally, if you have time1 , you might consider making the electric field explicitly time dependent(1.54)and investigate the behaviour of the system as a function of frequency.1.9.4 Your ReportOIn your report you should describe how you have chosen the algorithm and the value of and show thatthe method you have chosen is stable for this sort of equation and for the range of parameters you havechosen.You should think about the physics of the problem and try to identify the physically interesting parameter combinations. In particular you should ask yourself whether there is a qualitative change of behaviourat some value of the parameters or whether there is a value at which some special behaviour might beobserved.
Let your physical intuition guide your choice of parameters. If you can’t identify the physicallyinteresting values, then start by doing a broad sweep of the meaningful parameters to try to identify anyinteresting features. Once you have found the interesting parameter ranges you can look at them moreclosely.Try not to make meaningless changes. If the results can only depend on a particular combination ofparameters it is pointless to vary them individually.Your report should contain a representative selection of results and a discussion of the physics whichthey illustrate.Please include a listing of your program as an appendix well as sending it as an e-mail attachment toComputational-Physics, A. MacKinnon or ph.compphys@ic.ac.uk.
Do not send theseto my personal e-mail address.1 Thereare no marks for this sectionChapter 2Partial Differential Equations2.1 Types of EquationsThe PDE’s which occur in physics are mostly second order1. The work in this section is also considered inchapter III of Potter (1973) and chapter 17 of Press et al. (1992).For linear equations in 2 dimensions there is a simple classification in terms of the general equation0½0Y½½½0½ r K 0 f r K M r 0 3 ¾=r K 3 ¿£r j ½ x, rr r rrr(2.1)as shown in the following table 0À0 MM 0Ál M ConditionTypeEllipticHyperbolicParabolicExampleLaplace’s equation (2.2)Wave equation (2.3)Diffusion/Schrödinger equation (2.4)These are listed in their simplest form as follows (with the substitution0½ 0½r K 0 Ãr 00r ½ $ r 0 ½r K 0 M0 r 0r 0½ $ r½Ä rK0 r rrLaplace’s equationWave equationDiffusion equation ÂF where appropriate)(2.2)(2.3)(2.4)We shall consider each of these cases separately as different methods are required for each.2.2 Elliptic Equations — Laplace’s equationWe shall deal with elliptic equations later when we come to consider matrix methods (section 3).
For themoment it suffices to note that, apart from the formal distinction, there is a very practical distinction to bemade between elliptic equations on the one hand and hyperbolic and parabolic on the other hand. Generally speaking elliptic equations have boundary conditions which are specified around a closed boundary.Usually all the derivatives are with respect to spatial variables, such as in Laplace’s or Poisson’s Equation.Hyperbolic and Parabolic equations, by contrast, have at least one open boundary. The boundary conditions for at least one variable, usually time, are specified at one end and the system is integrated indefinitely.1 Thereare some exceptions such as the 4th order equation which describes torsion of a beam.13Partial Differential Equations14Thus, the wave equation and the diffusion equation contain a time variable and there is a set of initial conditions at a particular time.These properties are, of course, related to the fact that an ellipse is a closed object, whereas hyperbolæand parabolæ are open.2.3 Hyperbolic Equations — Wave equationsThe classical example of a hyperbolic equation is the wave equation0J0r 0 $ M0 r K 0rrThe wave equation can be rewritten in the formÅr M r Ki Ir $r r ror as a system of 2 equations [M r K; , r(2.5) M r K $ M r r K [rNote that the first of these equations (2.3a) is independent ofr rr r(2.6)(2.7)(2.8)and can be solved on it’s own.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
