Thompson - Computing for Scientists and Engineers (523188), страница 55
Текст из файла (страница 55)
We make the calculations over thesame range of x as previously (0 to 4), but we anticipate that the Noumerov methodwill be much more accurate than both the Euler-type methods and the iteration-integration Adams methods for first-order equations. We therefore begin with twice aslarge a stepsize as used previously, that is h is increased from 0.05 to 0.1. Thetests with the exponential-function solutions are straightforward to make, so I suggest that you make them yourself.Exercise 8.26(a) Compile Numerical DE_2; Noumerov Method so that FUNC returnsvalue = 1 and ANALYT returns value = exp (x) .
Then run the programfor the Noumerov method with xmin = 0, xmax = 4, and stepsize h = 0.1.Input the initial values y [1] = y (xmin) = 1 and y [2] = y (xmin + h)= exp (0.1) = 1.105171. Examine the computed y [k] and prepare a graph tocompare with the curves in Figure 8.15, which were computed using Euler-typemethods. How much more accurate is the Noumerov method than the most-accurate (method 3) Euler method? (Take into account that you are using twice thestepsize used in the Euler methods.)(b) Repeat the above calculations for step sizes of h = 0.05 and 0.025, goingout to x = 4 each time.
After taking into account that the number of steps toreach a given x is inversely proportional to h, do your estimates of the error perstep agree in order of magnitude with the estimate derived in Section 8.3 for theNoumerov method? n298SECOND-ORDER DIFFERENTIAL EQUATIONSFIGURE 8.17 Noumerov-method solution of the second-order harmonic-oscillator differentialequation with cosine solution. The solid curve is the analytical solution, the dashed curve is thescaled error in the numerical solution for stepsize of 0.1, and the dotted curve is the scaled error forstepsize of 0.2.We now return to the cosine function, examined previously for first- and second-order differential equations by the Euler methods (Figures 7.7 and 8.16). Wemake calculations over the same range of x as previously (0 to 7), but we anticipatethat the Noumerov method will be much more accurate than both the Euler-typemethods and the iteration-integration Adams methods.
Therefore we begin withtwice as large a stepsize as used previously, that is, h is increased from 0.05 to 0.1.The resulting values of y [k] and the scaled errors from the numerical integration areshown for stepsizes of 0.1 and 0.2 in Figure 8.17.Note the improvement by about five orders of magnitude in the accuracy of thecomputation, for very little increase in complexity of the algorithm, from Eulermethods to Noumerov’s method. To be fair in the general comparison, we shouldpoint out that the Noumerov algorithm is specially adapted to linear equations withthe first derivative eliminated, whereas the Euler-type methods can also handle nonlinear equations.Just to show that all this improvement isn’t just smoke and mirrors, why don’tyou verify for yourself the results shown in Figure 8.17?Exercise 8.27(a) Rename then recompile prognrm Numerical DE_2; Noumerov Method sothat FUNC now returns value -1 and ANALYT returns value cos (x) .
Run themodified program with xmin = 0, xmax = 7, and stepsize h = 0.1. Input thevalues y (xmin) = 0, y' (xmin) = 0. Examine the values of the computedy [k] and the percentage relative errors in order to verify approximate agreementwith the curves in Figure 8.17.8.5 NOUMEROV METHOD FOR LINEAR SECOND-ORDER EQUATIONS299(b) Repeat the calculations for stepsize h = 0.2, going out to x = 7. Aftertaking into account that the number of steps to reach a given x is inversely proportional to h, show that your estimates of the error per step agree with the0 (h6) estimates derived for the Noumerov method in Section 8.3.
nNow that we have a very accurate and efficient method for solving linear secondorder differential equations, it is worthwhile to illustrate its use in quantum mechanics, where the Schrödinger equation is a linear second-order equation of wide applicability for describing molecular, atomic, and nuclear phenomena.The quantum harmonic oscillatorThe one-dimensional harmonic oscillator in quantum mechanics is of basic importance, analogously to the simple pendulum in classical mechanics, as discussed (forexample) in Chapter 5 of the quantum-mechanics text by Cohen-Tannoudji et al.
Itis therefore interesting to apply the Noumerov algorithm to solve the Schrödingerequation numerically in order to verify the mathematical derivations of the quantizedenergy levels that are presented in quantum mechanics textbooks. Here we summarize the relevant mathematical results that are used for comparison with our numericalresults. Two textbooks and software that are of interest for computer-based quantum mechanics are those by Brandt and Dahmen.The time-independent Schrödinger equation for the wave function (x) in onedimension for an oscillator (Hooke’s-law) potential is written in terms of the particlemass m, the position coordinate x, and the oscillator angular frequency as(8.81)where h is Planck’s constant and E is the total energy.
For time-independent statesthe energies are quantized according to(8.82)where n is a non-negative integer (0,l,...) andIt is most convenient torewrite this equation in terms of dimensionless variables, which also provides numbers that are conveniently scaled for computation,Exercise 8.28(a) Consider in the Schrodinger equation the change of variables(8.83)(the denominator istimes the “oscillator length”), the scaling of the energy300SECOND-ORDER DIFFERENTIAL EQUATIONS(8.84)and renaming of the dependent variable fromto y. Show that (8.81) becomes(8.85)where the function F that is appropriate for the Noumerov algorithm is given by(8.86)in which the variables x and E are dimension-free.(b) Given the rule (8.82) for the quantized energies, show that in (8.86) thescaled energy is given by(8.87)with n = O,l,... . nIn the quantum mechanics of bound systems the signature of a time-independentstate, also called an “energy eigenstate,” is that for large x the value of y (x) tends tozero.
Therefore the test that we will make is to use the Noumerov algorithm to solve(8.85) with (8.86) numerically for values of E near the En. Energy eigenstatesshould exhibit the claimed behavior for large x.Analytic expressions for y (x) can be obtained for the energy eigenstates. Appropriately scaled solutions of (8.85), yn(x), that are of order of magnitude unity,are given by(8.88)where the Hermite polynomials of order n, Hn(x), may be obtained from(8.89)(8.90)Since the solutions are polynomials times a Gaussian, they eventually die off forlarge x approximately as a Gaussian, with the falloff being slower for the larger nvalues associated with larger total energy.Now that we have the analytic solutions of the differential equation (8.85) at energy eigenstates to check the accuracy of the Noumerov solutions, we are ready toadapt the program Numerical DE_2; Noumerov Method to solve (8.85).8.5NOUMEROV METHOD FOR LINEAR SECOND-ORDER EQUATIONS301Noumerov solution of the quantum oscillatorIn order to solve the modified Schrödinger equation (8.85) for the quantum oscillator, our program for the Noumerov algorithm given two subsections above shouldbe modified slightly in several places.(1) First, change the output file name to NUM.SE or some distinguishing name.(2) The next modification is to set xmin = 0 rather than making it an input variable, since the differential equation solution will start at zero.
To avoid codingerrors, the variable xmin should remain in the program.(3) A convenient change is to enclose all the main program below the input of newvariable n , xmax, and h within a while loop controlled by an input variablecalled energy for the scaled value of the variable E in (8.85). The programshould stay in the loop as long as energy > 0. The program can then be runfor a range of energy values near the scaled energies E,, given by (8.87).(4) Starting values should be obtained from the analytical values, rather than input.This is not essential, but it establishes a scale for they values which cannot otherwise be determined because we are solving a linear equation.
Therefore, sety[l] = ANALYT(xmin,n); y[2] = ANALYT(xmin+h,n);(5) Modify the coding and reference of the function F(x), called FUNC in the program, as follows:double FUNC (x, energy)/* d2y/dx2 = FUNC(x,energy) */double x, energy;double value;/* Schroedinger equation second derivative */value = x*x-energy;return value;(6) Modify the references to ANALYT, the analytical solution of the Schrödingerequation for quantum number n, from ANALYT (x) to ANALYT (x,n) . Rewritefunction ANALYT as follows:double ANALYT(x,n)/* Analytical solution of Schroedinger equation */double x;int n; /* energy level; n = 0,1,2,...