Thompson - Computing for Scientists and Engineers (523188), страница 56
Текст из файла (страница 56)
*/doubleint k;value,Hzero,Hone,Gauss,norm,Hn;302SECOND-ORDER DIFFERENTIAL EQUATIONSHzero = 1; Hone = 2*x; Gauss = exp ( -x*x/2 ) ;if (n==O) return Hzero*Gauss;if (n==l) return Hone*Gauss/sqrt(2);norm = 2;for (k=2; k<=n; k++) /* Hn by recurrence on k */{Hn = 2*(x*Hone-(k-l)*Hzero);Hzero = Hone; Hone = Hn; /* update */norm = 2*k*norm;return}Hn*Gauss/sqrt(norm);The two lowest-order Hermite polynomials and the recurrence relations for thehigher-order polynomials are as given by (8.89) and (8.90). The factor normcomes from (8.88).Having made these changes, with appropriate declarations of new variables, weare ready to use the program to explore the numerical solutions of the quantum-oscillator wave functions.Exercise 8.29Run your modified program for the Noumerov method applied to the Schrödinger equation with input values n = 0, xmax = 3, h = 0.02, which are all dimensionless variables. Within the energy loop use values near the energyeigenvalue, E0 = 1.
Check that your results agree with Figure 8.18. nFIGURE 8.18 Numerical solution of the quantum harmonic oscillator Schrödinger equation forenergies near the n = 0 eigenstate, using the Noumerov algorithm with stepsize h = 0.02.8.5NOUMEROV METHOD FOR LINEAR SECOND-ORDER EQUATIONS303Note that unusually small values of the stepsize h are needed to accurately track thedecaying exponential or Gaussian function. Figure 8.18 displays the numericalwave function values y [k] that approximate y(x). The three values of the energy E are close to the energy eigenstate value of unity for n = 0. Note that (towithin ±O.002) the numerical solution that best fits the analytical solution differsfrom unity by 0.025. If you vary the stepsize h you will find that the discrepancyfrom unity decreases steadily as h decreases. Also, the decaying behavior of thenumerical solution is extremely sensitive to the relative accuracy of the two startingvalues.
Why not explore these two aspects of the numerical solution yourself?Exercise 8.30(a) Run your Schrödinger equation program as in Exercise 8.29, but now varyh from say 0.1 to 0.01 with (for minimum confusion) energy fixed at the energy eigenvalue of unity. Is the solution converging toward the analytical solution as h decreases?(b) Modify the program so that y[2] is 1 + 10-6 times the analytical value. Repeat the calculations made in Exercise 8.29. Verify that for large x this part-permillion change in a starting value is amplified about 1000 times when x = 3 isreached. nIt should now be clear to you that it is not good practice to use a numerical algorithm iteratively in a direction for which the function tends to decrease overall.
Forexample, the increasing exponential function investigated by the Noumerov algorithm in Exercise 8.26 has part-per-million accuracy, whereas the accuracy at xmaxfor a decreasing exponential, as investigated in the two preceding exercises, is onlypart-per-cent or worse. Therefore, iteration by numerical algorithms should, whenever feasible, be carried out in the direction of increasing magnitude of y.In our present example of the quantum oscillator one could begin with the Gaussian solutions at large x, then iterate inward.
Unless you had selected an energyeigenstate, the solution would diverge near the origin because of the presence of thatother solution of the second-order differential equation, the one irregular (divergent)near the origin. Such an iteration direction should give results that are much lesssensitive to the stepsize than is the iteration direction we are using. An alternativescheme for numerical solution of such a “stiff’ differential equation as we are tryingto solve is presented in Section 8.6.Finally, it is interesting to take a quantum leap and try the Noumerov method forthe quantum oscillator on a state of higher energy, for which there is an oscillatingregion of the wave function followed by a decaying region for larger x if E is nearlyan eigenstate value.
Try it and see.Exercise 8.31Run the Noumerov method program for the Schrödinger equation with inputvalues n = 3, xmax = 3, h = 0.02, which are all dimension-free variables.Within the energy loop use values near the energy eigenvalue, E3 = 7, according to (8.87). Your results should be similar to those in Figure 8.19. n304SECOND-ORDER DIFFERENTIAL EQUATIONSFIGURE 8.19 Numerical solution of the quantum harmonic oscillator Schrödinger equation forenergies near the n = 3 eigenstate, using the Noumerov algorithm with stepsize h = 0.02.In Figure 8.19 the value that best fits the analytical solution is at an energy = E of 7.06, close to the analytical value of 7.
The numerical solution is indistinguishable from the analytical solution within the resolution of the graph.Clearly the values of E = 6.70 or E = 7.40 do not give energy eigenstates.Now that you understand how to solve the Schrödinger equation numerically toreasonable accuracy by using the Noumerov method, you will probably find it interesting (depending on your level of experience with quantum machanics) to exploreother potentials both for bound states (as here) and for scattering states.The subject of numerical solution of ordinary differential equations is vast. Youcan learn more about the analytical and numerical aspects in, for example, Chapters 9 and 10 of Nakamura’s book and in Jain’s compendium of numerical methods.
Several programs in C are provided in Chapter 15 of the numerical recipesbook by Press et al.8.6INTRODUCTION TO STIFF DIFFERENTIAL EQUATIONSWhen applying the Noumerov algorithm to solve the differential equations for exponential-type functions in Section 8.5 we were struck by an unusual contrast.Namely, if we iterated the numerical solution in the direction that the function wasgenerally increasing, as in Exercise 8.29, then the algorithm was remarkably accurate, typically at least part-per-million accurate for a stepsize h = 0.1.
By comparison, if the solution was iterated so that the function was generally decreasing, as inExercises 8.30 and 8.31, then even for h = 0.02 only part-per-cent accuracy wasobtained.8.6INTRODUCTION TO STIFF DIFFERENTIAL EQUATIONS305Although we suggested at the end of Section 8.5 that one way around this problem is always to iterate in the direction of x such that the function is generally increasing in magnitude, this is often not possible or not practicable. Alternativeschemes are presented in this section. We first describe what is meant by a “stiff’differential equation, then we give two prescriptions that may remove the stiffness,due to Riccati and to Madelung.What is a stiff differential equation?Suppose that in the second-order differential equation(8.91)we have that F (x) >> 0 for a wide range of x.
For any small range of x values,small enough that F does not change appreciably over this range, the solutions of(8.91) will be of the form(8.92)in which A+ and A- are constants.Exercise 8.32(a) Verify the correctness of this approximate form of the solution of (8.91) bycalculating the second derivative of this expression with respect to x, includingthe variation of F with respect to x. Show that if(8.93)then the approximation is appropriate.(b) Show that this condition can be weakened to(8.94)by appropriately moving the x origin and adjusting A+ and A-. nIn many scientific and numerical applications the appropriate solution of (8.91)is the exponentially decaying one.
For example, we may have the boundary condition thatso that A+ = 0 is required. If inaccuracies of the numerical algorithm or of roundoff error allow a small amount of exponentially-increasing solution to insinuate itself into the decaying solution, this increasing solutionpart will usually quickly increase. The differential equation is then “stiff’ to solve,and its solution is an unstable problem in the sense discussed in Section 4.3.306SECOND-ORDER DIFFERENTIAL EQUATIONSOne way of appreciating this insidious behavior is to consider the effect of asmall error, on the right-hand side of (8.91). The error in y" (x) is then FSince F >> 0, this amplifies the curve strongly upward or downward, dependingThis problem was introduced in Exercise 8.30 (b)on whetherand the subsequent discussion.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
