c19-3 (779618)
Текст из файла
19.3 Initial Value Problems in Multidimensions853The correct way to difference Schrödinger’s equation [1,2] is to use Cayley’sform for the finite-difference representation of e−iHt , which is second-order accurateand unitary:e−iHt '(19.2.35)In other words,1 + 12 iH∆t ψjn+1 = 1 − 12 iH∆t ψjn(19.2.36)On replacing H by its finite-difference approximation in x, we have a complextridiagonal system to solve. The method is stable, unitary, and second-order accuratein space and time. In fact, it is simply the Crank-Nicholson method once again!CITED REFERENCES AND FURTHER READING:Ames, W.F.
1977, Numerical Methods for Partial Differential Equations, 2nd ed. (New York:Academic Press), Chapter 2.Goldberg, A., Schey, H.M., and Schwartz, J.L. 1967, American Journal of Physics, vol. 35,pp. 177–186. [1]Galbraith, I., Ching, Y.S., and Abraham, E. 1984, American Journal of Physics, vol. 52, pp. 60–68. [2]19.3 Initial Value Problems in MultidimensionsThe methods described in §19.1 and §19.2 for problems in 1 + 1 dimension(one space and one time dimension) can easily be generalized to N + 1 dimensions.However, the computing power necessary to solve the resulting equations is enormous. If you have solved a one-dimensional problem with 100 spatial grid points,solving the two-dimensional version with 100 × 100 mesh points requires at least100 times as much computing.
You generally have to be content with very modestspatial resolution in multidimensional problems.Indulge us in offering a bit of advice about the development and testing ofmultidimensional PDE codes: You should always first run your programs on verysmall grids, e.g., 8 × 8, even though the resulting accuracy is so poor as to beuseless. When your program is all debugged and demonstrably stable, then you canincrease the grid size to a reasonable one and start looking at the results.
We haveactually heard someone protest, “my program would be unstable for a crude grid,but I am sure the instability will go away on a larger grid.” That is nonsense of amost pernicious sort, evidencing total confusion between accuracy and stability. Infact, new instabilities sometimes do show up on larger grids; but old instabilitiesnever (in our experience) just go away.Forced to live with modest grid sizes, some people recommend going to higherorder methods in an attempt to improve accuracy.
This is very dangerous. Unless thesolution you are looking for is known to be smooth, and the high-order method youSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).1 − 12 iH∆t1 + 12 iH∆t854Chapter 19.Partial Differential EquationsLax Method for a Flux-Conservative EquationAs an example, we show how to generalize the Lax method (19.1.15) to twodimensions for the conservation equation∂u= −∇ · F = −∂t∂Fx ∂Fy+∂x∂y(19.3.1)Use a spatial grid withxj = x0 + j∆(19.3.2)yl = y0 + l∆We have chosen ∆x = ∆y ≡ ∆ for simplicity.
Then the Lax scheme isun+1j,l =1 n(u+ unj−1,l + unj,l+1 + unj,l−1 )4 j+1,l∆t nnnn(F−− Fj−1,l+ Fj,l+1− Fj,l−1)2∆ j+1,l(19.3.3)Note that as an abbreviated notation Fj+1 and Fj−1 refer to Fx, while Fl+1 andFl−1 refer to Fy .Let us carry out a stability analysis for the model advective equation (analogof 19.1.6) withFx = vx u,Fy = vy u(19.3.4)This requires an eigenmode with two dimensions in space, though still only a simpledependence on powers of ξ in time,unj,l = ξ n eikx j∆ eiky l∆(19.3.5)Substituting in equation (19.3.3), we findξ=1(cos kx∆ + cos ky ∆) − iαx sin kx∆ − iαy sin ky ∆2whereαx =vx ∆t,∆αy =vy ∆t∆(19.3.6)(19.3.7)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.
To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).are using is known to be extremely stable, we do not recommend anything higherthan second-order in time (for sets of first-order equations).
For spatial differencing,we recommend the order of the underlying PDEs, perhaps allowing second-orderspatial differencing for first-order-in-space PDEs. When you increase the order ofa differencing method to greater than the order of the original PDEs, you introducespurious solutions to the difference equations. This does not create a problem if theyall happen to decay exponentially; otherwise you are going to see all hell break loose!19.3 Initial Value Problems in Multidimensions855The expression for |ξ|2 can be manipulated into the form|ξ|2 = 1 − (sin2 kx∆ + sin2 ky ∆)1− (α2x + α2y )2(19.3.8)The last two terms are negative, and so the stability requirement |ξ|2 ≤ 1 becomes1− (α2x + α2y ) ≥ 02∆∆t ≤ √2(vx2 + vy2 )1/2or(19.3.9)(19.3.10)This is an example of the general result for the N -dimensional Courantcondition: If |v| is the maximum propagation velocity in the problem, then∆t ≤ √∆N |v|(19.3.11)is the Courant condition.Diffusion Equation in MultidimensionsLet us consider the two-dimensional diffusion equation,∂u=D∂t∂2u ∂2u+ 2∂x2∂y(19.3.12)An explicit method, such as FTCS, can be generalized from the one-dimensionalcase in the obvious way.
However, we have seen that diffusive problems are usuallybest treated implicitly. Suppose we try to implement the Crank-Nicholson schemein two dimensions. This would give usHere1 2 n+1n2 n2 n+12 nun+1j,l = uj,l + α δx uj,l + δx uj,l + δy uj,l + δy uj,l2α≡D∆t∆2∆ ≡ ∆x = ∆yδx2 unj,l ≡ unj+1,l − 2unj,l + unj−1,l(19.3.13)(19.3.14)(19.3.15)and similarly for δy2 unj,l . This is certainly a viable scheme; the problem arises insolving the coupled linear equations. Whereas in one space dimension the systemwas tridiagonal, that is no longer true, though the matrix is still very sparse. Onepossibility is to use a suitable sparse matrix technique (see §2.7 and §19.0).Another possibility, which we generally prefer, is a slightly different way ofgeneralizing the Crank-Nicholson algorithm.
It is still second-order accurate in timeand space, and unconditionally stable, but the equations are easier to solve thanSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).1− (cos kx∆ − cos ky ∆)2 − (αy sin kx ∆ − αx sin ky ∆)24856Chapter 19.Partial Differential Equations(19.3.13).
Called the alternating-direction implicit method (ADI), this embodies thepowerful concept of operator splitting or time splitting, about which we will saymore below. Here, the idea is to divide each timestep into two steps of size ∆t/2.In each substep, a different dimension is treated implicitly:un+1j,l1 n+1/2= unj,l + α δx2 uj,l+ δy2 unj,l21 n+1/2n+1/2= uj,l+ α δx2 uj,l+ δy2 un+1j,l2(19.3.16)The advantage of this method is that each substep requires only the solution of asimple tridiagonal system.Operator Splitting Methods GenerallyThe basic idea of operator splitting, which is also called time splitting or themethod of fractional steps, is this: Suppose you have an initial value equation ofthe form∂u= Lu∂t(19.3.17)where L is some operator. While L is not necessarily linear, suppose that it can atleast be written as a linear sum of m pieces, which act additively on u,Lu = L1 u + L2 u + · · · + Lm u(19.3.18)Finally, suppose that for each of the pieces, you already know a differencing schemefor updating the variable u from timestep n to timestep n + 1, valid if that pieceof the operator were the only one on the right-hand side.
We will write theseupdatings symbolically asun+1 = U1 (un , ∆t)un+1 = U2 (un , ∆t)···(19.3.19)un+1 = Um (un , ∆t)Now, one form of operator splitting would be to get from n to n + 1 by thefollowing sequence of updatings:un+(1/m) = U1 (un , ∆t)un+(2/m) = U2 (un+(1/m) , ∆t)···un+1 = Um (un+(m−1)/m , ∆t)(19.3.20)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.
Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).n+1/2uj,l19.4 Fourier and Cyclic Reduction Methods857For example, a combined advective-diffusion equation, such as∂u∂2 u∂u= −v+D 2∂t∂x∂x(19.3.21)un+1/m = U1 (un , ∆t/m)un+2/m = U2 (un+1/m , ∆t/m)(19.3.22)···un+1 = Um (un+(m−1)/m , ∆t/m)The timestep for each fractional step in (19.3.22) is now only 1/m of the full timestep,because each partial operation acts with all the terms of the original operator.Equation (19.3.22) is usually, though not always, stable as a differencing schemefor the operator L.
Характеристики
Тип файла PDF
PDF-формат наиболее широко используется для просмотра любого типа файлов на любом устройстве. В него можно сохранить документ, таблицы, презентацию, текст, чертежи, вычисления, графики и всё остальное, что можно показать на экране любого устройства. Именно его лучше всего использовать для печати.
Например, если Вам нужно распечатать чертёж из автокада, Вы сохраните чертёж на флешку, но будет ли автокад в пункте печати? А если будет, то нужная версия с нужными библиотеками? Именно для этого и нужен формат PDF - в нём точно будет показано верно вне зависимости от того, в какой программе создали PDF-файл и есть ли нужная программа для его просмотра.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
