Thompson - Computing for Scientists and Engineers (523188), страница 53
Текст из файла (страница 53)
The way todo this is indicated in the following exercise.Exercise 8.21(a) Expand each of yk+1 and yk-1 in Taylor expansions about yk in order toshow that choosing D = -1 in (8.75) eliminates all odd-order derivatives from it.Now add these two values and corresponding formulas for their second derivatives to obtain(8.76)and the recurrence for derivatives(8.77)(b) In (8.75), with the substitutions from the two preceding equations, equatecorresponding derivatives on the left- and right-hand sides to show that B = 2,A = -1/12, C = 5/6, with remainder term R - h 6 y ( 6 ) k / 2 4 0 .
nThe original differential equation can now be written(8.78)which holds for any ordinary differential equation. It becomes useful, in addition tobeing true, when we use the linearity condition on its right-hand side, as you mayverify as follows.Exercise 8.22Define the iteration quantity(8.79)and make a binomial expansion approximation for yk in the terms of (8.78) thatare of order h2 and higher. Thus show that (8.78) can be written as(8.80)where the remainder R' is of order h6. n8.4 PROGRAMMING SECOND-ORDER EULER METHODS287Formula (8.80), used with (8.79) for yk and with neglect of the remainder, iscalled Noumerov’s algorithm for advancing the solution of the differential equationfrom k to k + 1.
In Section 8.5 we develop a program containing this algorithm,then we test its accuracy and efficiency before applying it to numerical calculationsof quantum oscillator wave functions.8.4 PROJECT 8A: PROGRAMMING SECOND-ORDEREULER METHODSOur goal in this project is to convert the three Euler-type algorithms, derived in Section 8.3 and summarized in Table 8.3, into working programs, then to test andcompare the accuracy of these algorithms in solving two simple differential equations, those for the exponential and cosine functions.Programming the Euler algorithmsThe Euler-type algorithms are straightforward to program, provided that care istaken to code the formulas in the order indicated for each method.
Recall that this isimportant because the formulas are approximations for advancing the solution of theequation from x to x + h, rather than equations that determine the y values at thesetwo values of the independent variable. In this sense, the = sign in the algorithmsis like the assignment operator in programming languages, the = in C and Fortranand the := in Pascal.For flexibility in testing and using the program, the second-derivative function,F (x,y), is a separate function referred to by the main program that contains thethree Euler-type algorithms as choices.
The function name is FUNC. Similarly,the analytical solution to the second-order differential equation that has the same twoboundary conditions that you use in your numerical solution is to be programmed infunction ANALYT.The overall structure of the program Numerical DE_2; Euler-Type Methodsis quite conventional. It begins with the usual opening of a file, NUMDE_2xx, to bewritten over or added to in order to save the output for graphic display.
The suffix“xx” in the file name is to be substituted by the function being calculated. For usewith the exponential function the program listing has xx = exp. For the cosinefunction an appropriate suffix is xx = cos. In the program supplied you have tomake new versions of the program for each new suffix. You should probably modify the file handling appropriately for your computing environment to avoid thisawkwardness.The three Euler-type algorithms appear adjacently rather than being hidden intheir own functions. By this means, which is not generally good programmingstyle, you easily see the similarities and differences between the coding of the threealgorithms. If you wish to use them independently as functions, you will need carein moving the appropriate variables into the functions.288SECOND-ORDER DIFFERENTIAL EQUATIONSAfter the algorithm loop of the code comes the comparison with analytical values.
Keeping the differential-equation solver loop separate from the comparisonloop makes the program structure clearer. Notice that we calculate both actual errors(err) and percentage relative errors (relerr). The former are more appropriate forfunctions that oscillate in sign but are bounded, such as the cosine and sine functions.
The second method of showing errors is better for functions with a largerange of variation, such as the exponential function computed in the next subsection.The complete program for solving second-order differential equations by theEuler-type methods summarized in Table 8.3 is as follows. Note that FUNC andANALYT are appropriate for exponential functions, being coded as F(x,y) = y andy (x) = exp (x), respectively.PROGRAM 8.1 Euler-type algorithms for numerical solution of second-order differential equations.
The three methods are summarized in Table 8.3.#include <stdio.h>#include <math.h>#define MAX 201main(){Euler-Type Methods/* Numerical DE_2;/* for Second-Order Differential Equations */FILE *fout;FILE *fopen();doubley[MAX],yp[MAX],err[MAX],relerr[MAX];double xmin,xmax,h,x,Funkh;doubleFunx,Funxl2,dYkp1,Ykp1,analytical;int choice,Nk,k;char wa;double FUNC() ,ANALYT();*/printf("Numerical DE; Second-Order Euler-type Methods\n");printf("Write over output (w) or Add on (a): ");fout = fopen("NUMDE_2Eexp",&wa);scarlf("%s",&wa);choice = 1;while ( choice > 0 )printf("\nChoose Euler method (zero to quit):\n");printf("1, 2, or 3; ");scanf("%i",&choice) ;if ( choice == 0 )printf("\nEnd Numerical DE_2; Euler-Type Methods"); exit(O);8.4PROGRAMMING SECOND-ORDER EULER METHODSif ( choice < 0 || choice > 3 ) /* || is 'or' */{printf("!! choice=%i only 1,2,3\n",choice); exit(l);}printf("Input xmin, xmax, h\n");scanf("%le%le%le",&xmin,&xmax,&h) ;Nk = (xmax-xmin)/h+l.l;if(Nk>MAX-1){printf("!! # of steps, %i, > %i\n",Nk,MAX-1); exit(l);printf("Input y(xmin), y' (xmin):\n");scanf("%le%le",&y[1],&yp[1] );Euler-type algorithm; Iteration in x *//*x = xmin;for ( k = 1; k <= Nk-1; k++ ){Funkh = FUNC(x,y,k)*h;yp[k+l] = yp[k]+Funkh;switch (choice)case 1: /* First method */yp[k+l] = yp[k]+Funkh;y[k+l] = y[k] +yp[k]*h; break;case 2: /* Second method */yp[k+1] = yp[k]+Funkh;y[k+l] = y[k]+(yp[k]+yp[k+1])*h/2; break;case 3: /* Third method */y[k+1] = y[k]+yp[k]*h+Funkh*h/2;yp[k+1] = yp[k]+(Funkh+FUNC(x+h,y,k+1)*h)/2;break;x = x+h;}/* end k for loop & Euler-type algorithms */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]);289290SECOND-ORDER DIFFERENTIAL EQUATIONSfprintf(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 */} /* end choice while loop */double FUNC(x,y,k)/* d2y/dx2 = FUNC(x,y) at y = y[k] */double Y[];double x;int k;double value;/* Using working example of F = y for FUNC */value = y[k];return value;}double ANALYT(x)/* Analytical solution of the differential equation */double x;double value;/* Using working example of y = exp(x) for ANALYT */value = exp(x);return value;}The following exercise may be used to check the correctness of coding the program before you start exploring use of Numerical DE_2; Euler-Type Methodsfor solving differential equations, The checking also illustrates a three-stage strategyfor developing numerical algorithms and associated programs.
The three stages are:(1) Verify correctness of coding the algorithm by comparison against test exampleswhere it should be exact.(2) Check appropriateness and limitations of the coded algorithm against knownfunctions.(3) Use the program for solving analytically intractable problems numerically, guided by its limitations that have been investigated at stage (2).Coding verification, stage (l), is illustrated in the following exercise.8.4 PROGRAMMING SECOND-ORDER EULER METHODS291Exercise 8.23(a) Show that each of the three Euler-type algorithms for second-order differential equations is exact if the function F (x,y) is identically zero. Similarly,show that the second and third algorithms are exact and identical if F (x,y) is aconstant. If you have physics background, you can use the correspondence tothe equations of kinematics with constant acceleration (y is displacement, x istime, and F is acceleration) to establish these results.(b) Modify the function FUNC so that it returns zero, and modify the functionANALYT so that it returns x.
Run the program with xmin = 0 and initial conditions y [l] = y (0) = 0, yp [1] = y' (0) = 1. Verify that all three methods give results that are exact within computer roundoff errors, with the firstderivative having the constant value of unity.(c) Modify function FUNC so that it returns the constant value of 2, and modify function ANALYT so that it returns x2. Run the program with xmin = 0and initial conditions y [l] = y (0) = 0, yp [l] = y' (0) = 0. Verify thatchoice = 2 and choice = 3 methods give results that are exact withinroundoff errors.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
