MacKinnon - Computational Physics, страница 5

The density and current density could be of mass,wherecharge, energy or something else which is conserved. Here we shall use the words charge and current forconvenience. We consider here the 1D form for simplicity. The equation is derived by considering spaceas divided into sections of length . The change in the total charge in a section is equal to the total currentcoming in (going out) through its endsOKª> Þ XZ KÅL$ > š ÑßYr ªÞ Üj:r(2.34)It is therefore useful to re–express the differential equation in terms of the total charge in the section andthe total current coming in through each face, so that we obtain a discrete equation of the form€ € € €€€r à I€ e€ T^á $ I á S#T(2.35)rá S#T is the current through the boundary between partswhere à represents the total charge in part â and Iâ and â .This takes care of the spatial part.

What about the time derivative? We can express the physics thus:The change in the charge in a cube is equal to the total charge which enters through itsfaces.2.6.2 The Diffusion Equation€€€€á S#T $ŸãÅ" à S#T $ à &£[€Q to€Q the diffusion€Q equation in discrete formSubstituting this into the equation€ of continuitydirectlyQS=T €Q leadsàà O ]ã Ç à S=T $ f à à egTJÈ(2.37)In many cases the current through a face is proportional to the difference in density (or total charge) betweenneighbouring cubesI(2.36)Partial Differential Equations19€ €€ †ä à € €ä à #S Tw ä àwhich is of course our simple method (2.27) of solution.To check whether this algorithm obeys the conservation law we sum over all , asshould beconserved.

Note that it helps to consider the whole process as taking place on a circle as this avoidsproblems associated with currents across the boundaries. In this case (e.g.)and it iseasy to see that the conservation law is obeyed for (2.37).2.6.3 Maxwell’s EquationsLet us now try to apply the same ideas to Maxwell’s equations. In free space we haveÝ ¥ E $må r HrÝ ¥ H | r E [r(2.38)(2.39)In order to reverse the derivation of these equations we consider space to be divided into cubes as before.For (2.6.3a) we integrate over a face of the cube and apply Stokes’ theoremæ1çÝ ¥ E šYè S > E †š è l¯$ r æ ç å H šYè Sr(2.40)(2.41)and the integral of the electric field, E, round the edges of the face is equal to minus the rate of change ofthe magnetic flux through the face, i.e.

Faraday’s law. Here we can associate the electric field, E with theedges and the magnetic field with the face of the cube.In the case of the diffusion equation we had to think in terms of the total charge in a cube instead of thedensity. Now we replace the electric field with the integral of the field along a line and the magnetic fieldwith the flux through a face.Note also that we can analyse (2.6.3b) in a similar way to obtain a representation of Ampère’s law.2.7 Dispersion"<K : %&; Y¿ Ì ‚+é ? e Í ª ƒÆ M Ë . Let us substitute (2.42) intoÆ the Lax algorithmwhere *e ƒ Te ƒ²ê Se6 ¿YÌ ‚+é ? WYXZ Í ªÏ 0 6 ¿YÌ ‚‘é ? W Í ª Ï Ç ¿ ́Í^ë ª ¿ Ì ÍJë ª È $Let us return to the Lax method (section 2.3.2) for hyperbolic equations. The solution of the differentialequation has the form(2.42)M O S ª $ e ª È1ì [O K Ç ¿ ́Í^ë ¿ ÌÎÍ^ë(2.43)é¿Ì ë ? ³ ´nµ " Ë O KÑ&i$ ) M O K µ )E¹ " Ë O KÑ&^[(2.44)OµFrom this we can derive a dispersion relation, *»í Ë , for the discrete equation and compare the resultwiththe original differential equation.

Since, in general, * could be complex we write it as* î * }) M Ë andforcomparereal and imaginary parts on both sides of the equation to obtain¢³ ´nµ î O ¿ e1ï ë ? ³ ´nµ Ë O K(2.45)1eïMµ )ð¹ î O ¿ ë ? ¯$ O K µ )E¹ñË O K[(2.46)OBy cancelling the common factors we now obtainPartial Differential Equations20M O ò\ó " KÑ&OK ¹ ËO³ ´nµ 0 " Ë O KÑ& ŒTaking the ratio of these or the sum of their squares respectively leads to the equationsò\ó ¹ ".î O %&;e ï¿ 0 ë? M O 0 µ 0 " Ky&\[(2.48)O Kw )E¹ Ë OTheus that in general the phase velocity is not M , although for long wavelengths,KË O õôfirst ofandtheseî O Pequationsô , we tellsrecover the correct dispersion relationship. This is similar to the situation when(2.47)we compare lattice vibrations with classical elastic waves: the long wavelength sound waves are OK butthe shorter wavelengths deviate.and,,The second equation (2.7b) describes the damping of the modes.

Again for smallbut (e.g.) short wavelengths,, are strongly damped. This may be a desirable property as shortwavelength oscillations may be spurious. After all we should have chosen to be small compared with anyexpected features. Nevertheless with this particular algorithmis not damped. Other algorithms,such as Lax–Wendroff, have been designed specifically to damp anything with a.ö ×OKKö f O K OËOKËOK Ž – ¢O ¢T2.8 Problems1. Explain the difference between hyperbolic, parabolic and elliptic partial differential equations, andgive an example of each.

What is the important physical distinction between hyperbolic and parabolicequations, on the one hand, and elliptic equations on the other?2. Describe the von Neumann procedure for analysing the stability of partial differential equations.3. Describe the physical principle behind the Courant–Friedrichs–Lewy condition as applied to thenumerical solution of partial differential equations.0j f 0 j $ 0j r 0 r K 0 r Ka)r j r0 j r 0 jr 0 jb)÷ r Kø r 0 $ r K $ f r K 0 r 0j r 0 j ˆ r r 0j r j ªrK 0 r 0 r KN r ¿c)0r Yj 0r j $ f r 0Yj r jr f j r 0 r K 0 r Kž r Kø r d)r j r j r $}r j r f r j $ j µ "<K ‡ &r Œir r K  ar K Œ r r KU E) ¹ e)r r rr r r5. The equationj‡jr 5ã r K ‡r rcan be represented by the difference equationj Q ‚ðù S=T%ƒ j Q ‚+ù ƒ å ° j#Q‚+ù S 0ƒ $ j#Q‚+ù S#Tƒ j Q ‚ðù ƒ $Gj#Q‚ðù eƒT ± žå ã K O ‡ [ODerive the truncation error of this difference equation.

KWrite down an alternative difference equation which is 2nd order accurate in both O and O .4. Are the following equations hyperbolic, elliptic or parabolic?6. The Dufort–Frankel (2.30) scheme is a method for the solution of the diffusion equation.Show that the method is unconditionally stable.Discuss the advantages and disadvantages of this method.Lagrangian Fluid21Ä varies in space " Ä Ä "<Ky&„& is0jjÄ r K 0 r Ä K r K [r r rQS#TU$ Ä QeT j Q‚ðù S#ƒT $}j Q‚+ù egƒTf O K »É f O K ú7.

The diffusion equation in a medium where the diffusion constantjjr r K Œ Ä r KU orr r rShow that the difference equationj Q ‚ðù S=T%ƒ $Gj Q ‚ðù ƒ Q j Q‚+ù S#ƒT $ f j Q ‚+ù ƒ j Q‚+ù egƒTÄOOK0j£K is not conserved.is not conservative, i.e. ûjr r¼Œ ÄConstruct an alternative difference scheme which is conservative.8. Show that the Lax (2.21) scheme for the solution of the advection equation is equivalent toj L $rrjKÑ0 $ T 0 0j[OÆ rr Kø Œ f O 0 Æ O rr K 0 higher order termsjüJg! " ) " Ë K$ * %&„& in the Lax scheme and explainExamine the behaviour of wave–like solutionsthe behaviour in terms of diffusion.9.

Describe what is meant by numerical dispersion.10. Lax–Wendroff method consists of 2 steps, just like Runge–Kutta or Predictor–Corrector. It is givenby€Q €QTÈ $ M O S=T $ È(2.49)Æ €QS#TR’ 0 0 O K €QÇ S#Æ TR’ 0 Æ° Æ S#TR’ 0 $ Æ eTR’ 0 ±(2.50)"<K %&Draw a diagram to illustrate the way this algorithm operates on an : grid.Show that the algorithm is stable provided the Courant–Friedrichs–Lewy condition is obeyed.

fOK.Show that the algorithm tends to dampen waves with wavelength ö€QS#TR’S#€ TR’ 0 0Æ QS=TƀQT 0 € Ç S#T Æ Q $ M O KÆ O€Q2.9 Project — Lagrangian Fluid CodeM Âý7¢yþ Üwithout changing shape (so that a sine waveSmall amplitude sound waves travel at speedremains a sine wave). However, when the amplitude is large, the sound speed differs between the peaksand the troughs of the wave and the wave changes shape as it propagates.

This project is to simulate thiseffect in one-dimension. A difference method is used in which the system is split up into cells and thecell boundaries are allowed to move with the local fluid velocity (a so–called Lagrangian (section—2.4)method). You should not require any library routines to complete this project.2.9.1 The Difference EquationsThe difference equations are,where,Æ Q‚+ù S S ÿZ ÿZ ƒ Æ Q‚+ù S e ÿZ ÿZ ƒ " þ Q ‚ðù ƒ $ þ Q‚+ùS=Tƒ & O „7 /KQ‚+ù S S#ÿ T„ƒ KgQ‚+ù S ƒ ÿ Q‚ðù S S ÿ ÿZ ƒ O :ZZ Æ ZQ ‚+ù ƒ constant ¥ " Q‚+ù ƒ & ï [þÜ(2.51)(2.52)(2.53)Lagrangian Fluidand,K Q‚ðù S ƒ ÿ22Q‚+ù ƒ / 7 "]K Q‚+ù S ƒ ÿZ V$ K Q‚+ù e ƒ ÿZ & :ÜQ‚+ù S ƒ ÿ(2.54)ƒQ‚+ùƒQ ‚+ù are the densitiesandare the positions and velocities of the cell boundaries.andZZÆÜþand pressures in the cells. / is the mass in each cell.It is useful for the purpose of programming to redefine things to get rid of the various half integerindices.

Hence we can write the equations asQB ‚+ù S##S T%T„ƒƒ QB ‚+ù ƒƒ " Q ‚+ù S#ƒ T%$ƒ Q‚+ùS#Tƒ & O „7 /ÆKÑQB ‚+ù Æ1K QB ‚+ù þ QB ‚+ù O þQ‚+ù ƒ / 7"]K BQ Ƃ+ù ƒ $VK BQ‚+eù T ƒ & :Ü(2.55)(2.56)(2.57)which maps more easily onto the arrays in a program.2.9.2 Boundary ConditionsQS#TKÑQeT8S#TK1Use periodic boundary conditions to allow simulation of an infinite wave train. Note, though, the presenceandin (2.9.1).