c19-2 (779617), страница 2
Текст из файла (страница 2)
The replacement (19.2.22) with n → n + 1 leavesus with a nasty set of coupled nonlinear equations to solve at each timestep. Oftenthere is an easier way: If the form of D(u) allows us to integratedz = D(u)du(19.2.23)analytically for z(u), then the right-hand side of (19.2.1) becomes ∂ 2 z/∂x2 , whichwe difference implicitly asn+1n+1− 2zjn+1 + zj−1zj+1(∆x)2(19.2.24)Now linearize each term on the right-hand side of equation (19.2.24), for examplezjn+1 ≡ z(un+1) = z(unj ) + (un+1− unj )jj∂z ∂u j,n(19.2.25)= z(unj ) + (un+1− unj )D(unj )jThis reduces the problem to tridiagonal form again and in practice usually retainsthe stability advantages of fully implicit differencing.Schrödinger EquationSometimes the physical problem being solved imposes constraints on thedifferencing scheme that we have not yet taken into account.
For example, considerthe time-dependent Schrödinger equation of quantum mechanics. This is basically aparabolic equation for the evolution of a complex quantity ψ. For the scattering of awavepacket by a one-dimensional potential V (x), the equation has the formi∂2 ψ∂ψ= − 2 + V (x)ψ∂t∂x(19.2.26)(Here we have chosen units so that Planck’s constant h̄ = 1 and the particle massm = 1/2.) One is given the initial wavepacket, ψ(x, t = 0), together with boundarySample 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).The Crank-Nicholson method can be generalized similarly.The second complication one can consider is a nonlinear diffusion problem,for example where D = D(u).
Explicit schemes can be generalized in the obviousway. For example, in equation (19.2.19) write852Chapter 19.Partial Differential Equationsconditions that ψ → 0 at x → ±∞. Suppose we content ourselves with firstorder accuracy in time, but want to use an implicit scheme, for stability.
A slightgeneralization of (19.2.8) leads to"#"#n+1n+1ψj+1− 2ψjn+1 + ψj−1=−+ Vj ψjn+1(∆x)2(19.2.27)for whichξ=1k∆x4∆t2+ Vj ∆t1+isin(∆x)22(19.2.28)This is unconditionally stable, but unfortunately is not unitary. The underlyingphysical problem requires that the total probability of finding the particle somewhereremains unity. This is represented formally by the modulus-square norm of ψremaining unity:Z∞−∞|ψ|2 dx = 1(19.2.29)The initial wave function ψ(x, 0) is normalized to satisfy (19.2.29). The Schrödingerequation (19.2.26) then guarantees that this condition is satisfied at all later times.Let us write equation (19.2.26) in the formi∂ψ= Hψ∂t(19.2.30)where the operator H isH =−∂2+ V (x)∂x2(19.2.31)The formal solution of equation (19.2.30) isψ(x, t) = e−iHt ψ(x, 0)(19.2.32)where the exponential of the operator is defined by its power series expansion.The unstable explicit FTCS scheme approximates (19.2.32) asψjn+1 = (1 − iH∆t)ψjn(19.2.33)where H is represented by a centered finite-difference approximation in x.
Thestable implicit scheme (19.2.27) is, by contrast,ψjn+1 = (1 + iH∆t)−1 ψjn(19.2.34)These are both first-order accurate in time, as can be seen by expanding equation(19.2.32). However, neither operator in (19.2.33) or (19.2.34) is unitary.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).ψjn+1 − ψjni∆t19.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∆t.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
