Thompson - Computing for Scientists and Engineers (523188), страница 54
Текст из файла (страница 54)
Also verify that the fist derivative has the value 2x. nWith this program verification completed, we are ready for stage (2). That is,we investigate how well these algorithms handle typical differential equations whosesolutions are already known.Euler algorithms and the exponential functionThe second-order differential equation d2y /dx2 = y (x), with solution the increasing exponential function, y(x) = exp (x) if y (0) = y' (0) = 1, is a fairlydemanding differential equation against which to test the three Euler-type algorithmsin Table 8.3.
The same solution arises for the first-order differential equationdy/dx = y (x) with boundary condition y (0) = 1, as is investigated in Section 7.4, so a comparison of methods can be made. Indeed, method 1 for secondorder equations is identical to the forward-predictor method for first-order equations,because second-derivative terms in advancing y (x) are ignored in both algorithms.The three algorithms for the second-order equations are compared in Figure 8.15 for a step size h = 0.05. Method 1 produces an error of about 10% byx = 4, that is after 4/0.05 = 80 steps, or an error per step of about 1.25 × 10-3.Notice that the line labeled %err(l) in this figure is identical with that labeled%err(f) in Figure 7.5 because they involve the same numerical methods for handling the same differential equation.
Notice that the errors decrease by nearly an orderof magnitude as we go from method 1 to method 2 to method 3. The numerical effort involved hardly increases, as you can see by comparing the coding for the threechoices in Program 8.1, Numerical DE_2; Euler-Type Methods, which islisted in the preceding subsection.292SECOND-ORDER DIFFERENTIAL EQUATIONSFIGURE 8.15 Errors in three Euler-type methods for the second-order exponential equation.In this example of the increasing exponential function, all the derivatives areequal to the function value, so that once the function starts heading upward it continues uniformly upward. Even in method 3 derivatives beyond the second are ignored, so the error also increases uniformly, albeit slowly.
Because, for all three methods, we see from Figure 8.15 that the errors increase proportionally to x (since theerror lines are straight), the fractional error per step is constant. This result can beproved analytically for the exponential function. It is also of interest to investigatenumerically the dependence of the errors in each method on the stepsize h.Exercise 8.24(a) Compile Numerical DE_2; Euler-Type Methods so that FUNC returnsvalue y [k] and ANALYT returns value exp (x) .
Then run the program foreach of the three Euler methods with xmin = 0, xmax = 4, and stepsizeh = 0.05. Input y (xmin) = y' (xmin) = 1. Examine the values of the computed y [k] and the percentage relative errors in order to verify approximateagreement with the curves in Figure 8.15.(b) Repeat the above calculations for step sizes of h = 0 .1 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 with the estimates summarized in Table 8.3?This exercise should convince you that the Euler-type algorithms for second-orderdifferential equations are fairly inaccurate, because the error to propagate the solutionto a given x decreases one power of h more slowly than the error per step that is indicated in Table 8.3.
Probably only method 3 is suitable for typical applications.The merits of the three methods when used for oscillatory functions are explored inthe next subsection.8.4PROGRAMMING SECOND-ORDER EULER METHODS293FIGURE 8.16 Errors in three Euler-type methods for the second-order harmonic equation withsolution the cosine function.Euler algorithms and the cosine functionFor oscillatory functions it is interesting to examine the buildup and decay of the errors of the three Euler-type methods as the function varies. The function that exactlysolves the second-order linear differential equation d2y ldx2 = -y (x) with boundary conditions y (0) = 1 and y' (0) = 0, namely y(x) = cos (x), is the same asexplored for the first-order differential equation solvers at the end of Section 7.5.There the equation was dyldx = -sin (x) with y (0) = 1.
In real applications itwould not be sensible to go to the extra trouble of solving a second-order equation ifthere is a first-order equation available, but doing so here provides insight from thecomparison. We shall see immediately why there is extra trouble.The numerical solutions using the three Euler-type algorithms with a stepsize ofh = 0.05 are shown in Figure 8.16 and may be compared directly with the first-order solutions shown in Figure 7.7. Notice that only one of the second-order methods we use is more accurate than any of our first-order methods, namely method 3.It is more accurate than the first-order forward-difference method by about a factorof 100 and it is of comparable accuracy to the first-order central-difference and trapezoid methods.
It is worthwhile to convince yourself of the correctness of the numerical results for the second-order equation algorithms.Exercise 8.25(a) Rename and recompile program Numerical DE_2; Euler-Type Methodsso that FUNC now returns value - y [k] and ANALYT returns value cos (x) .Run the modified program for each of the three Euler methods with xmin = 0,xmax = 7, and a stepsize h = 0.05. Input the values y (xmin) = 1,y' (xmin) = 0. Examine the values of the computed y [k] and the percentageSECOND-ORDER DIFFERENTIAL EQUATIONS294relative errors in order to verify approximate agreement with the curves in Figure 8.16.(b) Repeat these calculations for stepsizes of h = 0.1 and 0.025, going out tox = 7 each time.
(The array sizes will need to be at least MAX = 282 for thelatter calculation.) After taking into account that the number of steps to reach agiven x is inversely proportional to h, do your estimates of the error per stepagree with the estimates summarized in Table 8.3? nWe may conclude from these numerical studies that numerical solutions of low-orderdifferential equations are usually more accurate than the solutions obtained fromhigher-order equations.
The only advantage of using a second-order equation over afirst-order equation for the cosine function example is that one obtains numerical estimates of slopes as a by-product of estimating the values of the equation solution.In the Euler-type methods investigated in this section and in Sections 7.4 and7.5 there was no restriction to linear differential equations, although our examplesdo not illustrate this. If the restriction to linear equations can be made, so thatF (x,y) = F (x)y in (8.54), there remains a very large range of scientifically interesting differential equations, to be solved numerically.
The more accurate methods,such as Noumerov’s algorithm that we derived in Section 8.3, are then appropriate.The numerics of this algorithm is the topic of the next section.8.5PROJECT 8B: NOUMEROV METHOD FOR LINEARSECOND-ORDER EQUATIONSIn this project our goal is to convert the Noumerov algorithm derived in Section 8.3into a working program, to evaluate the method, then to apply it to a problem inquantum mechanics. We first outline the program structure and list the program,then provide two test functions (exponential and cosine) for evaluating the method.In the last two subsections we apply the method and adapt the program to solve numerically the quantum-mechanical Schrödinger equation for a harmonic oscillator inone dimension.Programming the Noumerov methodCoding and testing of the Noumerov algorithm, (8.79) and (8.80), is quite straightforward.
For ease of presentation it is coded in-line, but it is clearly distinguishedfrom the main body of the program. The function that is iterated in-place is Y in(8.79), whereas it is the desired solution, y, that is stored in an array. The startingof the algorithm is assumed to be done by providing the first two values of y, iny [1] and y [2] . In some uses it might be more convenient to have the initial valueof y and its first or second derivatives.
These can be interconverted by use of aTaylor expansion and of the differential equation.8.5NOUMEROV METHOD FOR LINEAR SECOND-ORDER EQUATIONS295The overall program structure is by now familiar to you. The program Numerical DE_2; Noumerov Method begins with a file-control section for preparing fileNUMDE_2 . N, which receives the output for graphical or other subsequent processing.The file name should probably be made different for each application of the program, for example, the exponential, cosine, or harmonic-oscillator functions investigated below. The range of x values and the stepsize are requested, together withtwo initial values of y.
Then the algorithm discussed above is executed, followed bycomparison with analytic values made in terms of actual errors (err) and percentagerelative errors (relerr).The function FUNC is the s-dependent function that returns the second derivative, F(x). The function ANALYT should be coded to return the analytical solutionof the differential equation that has the same boundary conditions that the startingvalues of y indicate. It will be least confusing if you make new versions of the program whenever you change the file name or either of these function names.
Theversionof Numerical DE_2; Noumerov Method that we now show is setup fortesting the evolution of the cosine solution.PROGRAM 8.2 Noumerov algorithm for second-order linear differential equations without firstderivative terms.#include <stdio.h>#include <math.h>#define MAX 201main()/* Numerical DE_2;Noumerov Method*//* for Second-Order Differential Equations */FILE *fout;FILE *fopen();doubley[MAX],err[MAX],relerr[MAX];double xmin,xmax,h,h2,Y1,Y2,dYk,Yk,x;doubleFunx,Funx12,dYkp1,Ykp1,analytical;int Nk,k;char wa;double FUNC(),ANALYT();printf("Numerica1 DE; Second Order: Noumerov Method\n");printf("Write over output (w) or Add on (a): ");scanf("%s",&wa);fout = fopen("NUMDE_2.N",&wa);printf("Input xmin, xmax, h\n");scanf("%le%le%le",&xmin,&xmax,&h);Nk = (xmax-xmin) /h+1.1;if(Nk > MAX-1)296SECOND-ORDER DIFFERENTIAL EQUATIONS{printf("!! # of steps, %i, > %i\n",Nk,MAX-1); exit(l);}printf("Input y(xmin), y(xmin+h): ");scanf("%le%le",&y[1],&y[2]);/*Noumerov algorithm; Iteration in x */h2 = h*h;Y1 = (1-h2*FUNC(xmin)/l2)*y[1]; /* starting Y values */Y2 = (1+h2*FUNC(xmin+h)/l2)*y[2];dYk = Y2-Y1;Yk = Y2;x = xmin+2*h;for ( k = 2; k <= Nk-1; k++ )Funx = h2*FUNC(x); Funxl2 = (l+Funx/l2);dYkp1 = dYk+Funx*y[k]; /* increment in Y */Ykp1 = Yk+dYkpl; /* new Y value */y[k+1] = Funxl2*Ykpl; /* next y value */x = x+h; dYk = dYkp1; Yk = Ykp1; /* recycle values */}/* end k for loop &end Noumerov algorithm */printf("\nComparison with analyticalx = xmin;for ( k = 1; k <= Nk; k++ )values\n");analytical = ANALYT(x); /* analytical and error */err[k] = analytical-y[k];relerr[k] = lOO*err[k]/analytical;printf("%6.2le %10.6lg %10.6lg %10.6lg\n",x,y[k] ,err[k] ,relerr[k]) ;fprintf(fout, "%6.2le %10.6lg %10.6lg %10.6lg\n",x,y[k] ,err[k] ,relerr[k]) ;x = x+h; /* next x */} /* end k output loop */printf("\n End of Numerical DE_2; Noumerov Method");double FUNC(x)/* d2y/dx2 = FUNC(x) *//* Noumerov method is second-order and linear in y */double x;{double value;/* Using working example of F = -1 for FUNC */value = -1;8.5NOUMEROV METHOD FOR LINEAR SECOND-ORDER EQUATIONS297return value;double ANALYT(x)/* Analytical solution of differential equation */double x;double value;/* Using working example of y = cos(x) for ANALYT */value = cos(x);return value;With the program coding completed, we are ready to test it for correctness andaccuracy in the next subsection, then apply it in the final subsection.Testing Noumerov for exponentials and cosinesWe test the accuracy and efficiency of the Noumerov algorithm for second-order linear differential equations by using the same test functions, exponentials and cosines,as in the previous tests in Chapters 7 and 8.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
