Thompson - Computing for Scientists and Engineers (523188), страница 47
Текст из файла (страница 47)
It would be desirable to follow the extrapolation (predictor)step by a consistent corrector step in which the original differential equation is usedto improve the solution accuracy. Such predictor-corrector methods are described inChapter 8 of Vandergraft and in the book by Zill. Jain’s treatise on numerical solution of differential equations contains much more than you will probably ever use tomeet your needs.In the next section we describe a program that you can use to explore numerically the four solution methods that we have developed.7.5PROJECT 7: PROGRAM FOR SOLVINGFIRST-ORDER EQUATIONSIn the preceding section we worked out predictor formulas for numerical solution offirst-order differential equations of the kind (7.34).
We also gained experience withthe Euler forward- and central-predictor methods, anticipating the program that ispresented in this section. We now supplement these methods with the numerics andcoding of the Adams integration formulas, (7.46) and (7.47), derived in the preceding subsection. We first describe the parts of the program that relate to the Adamsmethods, we describe how to test the methods, then we explore use of the programfor solving differential equations numerically.Programming the differential equation solverThe Adams predictor methods involve some extra steps in programming and codingthan are needed for the Euler predictors, so we now describe these steps.The realities of computing require us to be more specific about the iteration procedures (7.46) and (7.47), because we must have a stopping criterion. In the program Numerical DE_1 for cases 3 and 4 (the Adams predictors) only one of twostopping criteria needs to be satisfied for the iteration on n to quit. The first criterionis that no more than nmax iterations are made, where nmax is an input parameterthat forces a definite stopping of the iteration.
The alternative criterion is that thefractional change in the y value since the last iteration should be less than an inputvalue epsilon.Exercise 7.20(a) Explain why it is both inefficient and uninformative to use nmax as the onlystopping criterion in the iterations (7.46) and (7.47).(b) Explain why it is dangerous (for the computer) to use only epsilon as theiteration-stopping criterion. n248INTRODUCTION TO DIFFERENTIAL. EQUATIONSSince either criterion used alone is ineffective, but their combination is informativeand prudent, in Numerical DE_1 if either criterion is satisfied, the iteration stops.This explains in case 3 and 4 the while ( && ) statements for the iterationloops, where && is C-language for the logical “and” operation.The comparison with epsilon is straightforward for the Adams-trapezoid case(3), but for Adams-Simpson (case 4) it is more complicated because we are iteratingon two adjacent k values.
A reasonable compromise that is made in the program isto average the absolute fractional differences for these two values. Experienceshows that if only one difference is used and the fractional difference for the other isignored, then the error for the ignored k values tends to grow unreasonably quickly.Averaging forces a moderating compromise.In either Adams method, if epsilon is made comparable to the number of significant figures that your computer is using, there will be a complicated and confusing interplay between convergence of the iteration and computer limitations fromsubtractive cancellation and other roundoff problems of the kinds examined inSection 4.3.So, with all these warnings we should be ready to plunge into the maelstrom ofnumerical mathematics.
Here is the complete program.PROGRAM 7.2 Numerical solution of first-order differential equations by Euler and Adamspredictor methods.#include <stdio.h>#include <math.h>#define MAX 201/* Numerical DE_l; Euler & Adams Predictors */for First-Order Differential Equations *//*FILE *fout;FILE *fopen();doubley[MAX],err[MAX],relerr[MAX];double xmin,xmax,h,epsilon,x,ylast1,ylast2,absdiff,analytical;intchoice,Nk,nmax,k,kmin,kstep,n;char wa;double FUNC(),ANALYT();printf("Numerica1 DE; First Order\n");printf("Write over output (w) or Add on (a):\n");fout = fopen("NUMDE_1",&wa);scanf("%s",&wa);choice = 1;while ( choice > 0 )printf("\nChoose method of solution (zero to quit):\n");7.5PROGRAM FOR SOLVING FIRST-ORDER EQUATIONSprintf("1 Euler-forward predictor;\n");printf("2 Euler-central predictor;\n");printf("3 Adams-trapezoid predictor;\n");printf("4 Adams-Simpson predictor;\n");scanf("%i",&choice) ;if ( choice == 0 ){printf("\nEnd Numerical DE; First Order"); exit(O);}if ( choice < 0 || choice > 4 ){printf("!! choice=%i only 1,2,3 or 4\n",choice);exit(l);}printf("Input xmin, xmax, h\n");scanf ("%1e%1e%1e",&min,&xmax,&h) ;Nk = (xmax-xmin)/h+l.l;if(Nk>MAX-1){printf("!! # of xsteps, %i, > %i\n",Nk,MAX-1); exit(l);}/* Boundary values */switch (choice){case 1: /* Euler-forward predictor */{printf("1: Input y(xmin);\n");scanf("%le",&y[l]);kmin=l; kstep = 1; break;}case 2: /* Euler-central predictor */{printf ("2: Input y(xmin),y(xmin+h);\n");scanf ("%le%le",&y[l],&y[2]);kmin=2; kstep = 2; break;}case 3: /* Adams-trapezoid predictor */{printf("3: Input y(xmin),nmax,epsilon;\n");scanf("%le%i%le",&y[l],&nmax,&epsilon) ;kmin=l; kstep = 1; break;}case 4: /* Adam-Simpson predictor */{printf("4: Inputy(xmin),y(xmin+h),nmax,epsilon;\n");scanf("%le%le%i%le",&y[l],&y[2],&nmax,&epsilon);kmin=2; kstep = 2; break;249250INTRODUCTION TO DIFFERENTIAL EQUATIONS/* Recurrence in x */x=xmin;for ( k = kmin k < Nk; k = k+kstep ){switch (choice){case 1: /* Euler-forward predictor */{y[k+l] = y[k]+h*FUNC(x,y,k); break;}case 2: /* Euler-central predictor */{x = x+h;y[k+l] = y[k-1]+2*h*FUNC(x,y,k);y[k+2] = y[k]+2*h*FUNC(x+h,y,k+l); break;}case 3: /* Adams-trapezoid predictor */{ /* Start iteration */y[k+l] = y[k]; n = 1; absdiff = 1;while ( n <= nmax && absdiff > epsilon ) /* iterate */{ylast1 = y[k+1];y[k+1] = y[k]+(FUNC(x,y,k)+FUNC(x+h,y,k+l) )*h/2;absdiff = fabs( (y[k+1]-ylast1)/ylast1);n = n+l;}break;}case 4: /* Adams-Simpson predictor */x = x+h; /* Start iteration */y[k+1] = y[k]; y[k+2] = y[k]; n = 1; absdiff = 1;while ( n <= nmax && absdiff > epsilon ) /* iterate */ylast1 = y[k+l]; ylast2 = y[k+2];y[k+1] = y[k-1]+(FUNC(x-h,y,k-1)+4*FUNC(x,y,k)+FUNC(x+h,y,k+l))*h/3;y[k+2] = y[k]+(FUNC(x,y,k)+4*FUNC(x+h,y,k+l)+FUNC(x+2*h,y,k+2))*h/3;/* For convergence average fractional diferences */absdiff = (fabs((y[k+1]-ylast1)/ylast1)+fabs((y[k+2]-ylast2)/ylast2))/2;7.5PROGRAM FOR SOLVING FIRST-ORDER EQUATIONSn = n+l;}break;}}x = x+h;} /* end k for loop */printf("Comparison with analytical values\n");x=xmin;for ( k = 1; k <= Nk; k++ )analytical = ANALYT(x);/* analytical and error */err[k] = analytical-y[k];relerr[k] = lOO*err[k]/analytical;x = x+h;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]);} /* end k output loop */} /* end choice while loop */double FUNC(x,y,k)/* dy/dx = FUNC */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);251252INTRODUCTION TO DIFFERENTIAL EQUATIONSreturn value;}Exercise 7.21(a) Key in those parts of Numerical DE_-1 thatrefer to case 3 and case 4,assuming that you did Exercise 7.18.
(Otherwise, key in the complete program.) Note that the main program refers to f in (7.34) as FUNC . For use aftertesting you will substitute the coding for the function of interest to you. For ourworking example the function is y. Don’t anticipate the exact value by usingexp (x ) instead.(b) Run this version of the program for an interesting range of x values. Forsimplicity xmin = 0,xmax = 4.O and h = 0.05, as in Exercise 7.18, arereasonable choices. The file ( NUMDE_1 ) will record your output for use in aspreadsheet or graphics application.
The iteration-control values used in Figure7.6 are nmax = 5 and epsilon = 10-3, with the latter convergence criterion usually being satisfied after n = 3 iterations.(c) Compare the numerical values just computed with analytical values for thesolution of the differential equation that are computed by function ANALYT. RunNumerical DE_1 to check your output against the values shown graphically inFigure 7.6. nThe exponential solution of our working-example differential equation is shown inFigure 7.6 for the two Adams predictor methods, along with the relative errors forthe two integration methods (trapezoid and Simpson) shown as err(trap) anderr(Simp) respectively. Note the major error reduction, by about a factor of 100,for Adams-trapezoid method over Euler-forward method and for Adams-Simpsonover Euler-central. This is not all gain without pain, because if n iterations are usedin the Adams methods, the computing time is increased by about the same factor.Dependent on how complex it is to calculate FUNC f, it might sometimes be moreefficient for the same numerical accuracy to use a smaller stepsize, h, in one of theEuler methods.
Try it and compare.Exploring numerical first-order equationsNow that we have gone to the trouble to develop a practical program for numericalsolution of first-order differential equations, Numerical DE_1, it is interesting toexplore other equations with it.The first equation that I suggest trying is(7.48)This has the well-known exact solution(7.49)7.5PROGRAM FOR SOLVING FIRST-ORDER EQUATIONS253FIGURE 7.6 Adams-predictor solutions of the first-order differential equations (7.34), (7.35)with boundary condition y (0) = 1. The solid line shows the analytical solution, exp (x), thedashed line shows l02 times the percentage error in the Adams-trapezoid numerical solution, and the3dotted line shows 10 times the percentage error in the Adams-Simpson solution.Because of the oscillatory behavior and the initial decrease of the function as onemoves away from x = 0, we can see whether the error buildup is sensitive to theoscillatory behavior, something that we could not do for the exponentially increasingfunction used as the working example.For simplicity, starting at x = 0 and going to just over x =is a goodchoice.
This allows us to see how nearly periodic is the numerical solution of(7.48). Because the exact solution becomes zero if x is an odd multiple oftheoutput is less confusing if actual errors, ek from (7.32), are used instead of the relative errors, rk from (7.33), shown in the preceding computing exercises in this chapter.Exercise 7.22(a) In Numerical DE_1 modify function FUNC to assign -sin(x) to valuefor the first derivative, and also modify function ANALYT to assign cos (x) tovalue for the analytical solution of the differential equation (7.48).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
