Deturck, Wilf - Lecture Notes on Numerical Analysis (523142), страница 15
Текст из файла (страница 15)
The program# first calls eulermethod to obtain the array ynew of guessed values# of y at x+h. It then refines the guess repeatedly, using the trapezoidal# rule, until the previous guess, yguess, and the refined guess, ynew, agree# within a tolerance of eps in all components. Program then computes dist,# which is the largest deviation of any component of the final converged# solution from the initial Euler method guess.
If dist is too large66>>>>>>>>>>>>>>>>>The Numerical Solution of Differential Equations# the mesh size h should be decreased; if too small, h should be increased.ynew:=eulermethod(yin,x,h,f);yfirst:=ynew;ll:=nops(yin);allnear:=false;while(not allnear) doyguess:=ynew;ynew:=[];for i from 1 to ll doynew:=[op(ynew),yin[i]+(h/2)*(f(x,yin,i)+f(x+h,yguess,i))];od;allnear:=true;for i from 1 to ll do allnear:=allnear and abs(ynew[i]-yguess[i])<eps od:od; #end whiledist:=max(seq(abs(ynew[i]-yfirst[i]),i=1..ll));RETURN([dist,ynew]):end:The two programs above each operate at a single point x and seek to compute theunknowns at the next point x + h.
Now we need a global view, that is a program that willcall the above repeatedly and increment the value of x until the end of the desired rangeof x. The global routine also needs to check whether or not the mesh size h needs to bechanged at each point and to do so when necessary.>>>>>>>>>>>>>>>>>>>>>>>>>>trapglobal:=proc(f,y0,h0,xinit,xfinal,eps,nprint)local x,y,y1,h,j,arr,dst,cnt:# Finds solution of the ODE system y’=f(x,y), where y is an array# and f is array-valued. y0 is initial data array at x=xinit.# Halts when x>xfinal.
eps is convergence criterion for# trapezoidal rule; Prints every nprint-th value that is computed.x:=xinit:y:=y0:arr:=[[x,y[1]]]:h:=h0:cnt:=0:while x<=xfinal doy1:=traprule(y,x,h,eps,f):y:=y1[2]:dst:=y1[1];# Is dst too large? If so, halve the mesh size h and repeat.while dst>3*eps doh:=h/2; lprint(‘At x=‘,x,‘h was reduced to‘,h);y1:=traprule(y,x,h,eps,f):y:=y1[2]:dst:=y1[1];od:# Is dst too small? If so, double the mesh size h and repeat.while dst<.0001*eps doh:=2*h; lprint(‘At x=‘,x,‘h was increased to‘,h);y1:=traprule(y,x,h,eps,f):y:=y1[2]:dst:=y1[1];od:# Adjoin newly computed values to the output array arr.x:=x+h; arr:=[op(arr),[x,y[2]]]:# Decide if we should print this line of output or not.cnt:=cnt+1: if cnt mod nprint =0 or x>=xfinal then print(x,y) fi;2.11 Maple programs for the trapezoidal rule>> RETURN(arr);> end:67od:The above three programs comprise a general package that can numerically solve systemsof ordinary differential equations.
The applicability of the package is limited mainly by thefact the Euler’s method and the Trapezoidal Rule are fairly primitive approximations tothe truth, and therefore one should not expect dazzling accuracy when these routines areused over long intervals of integration.2.11.1Example: Computing the cosine functionWe will now give two examples of the operation of the above programs. First let’s computethe cosine function.
We will numerically integrate the equation y 00 + y = 0 with initialconditions y(0) = 1 and y 0 (0) = 0 over the range from x = 0 to x = 2π. To do this we usetwo unknown functions y1 , y2 which are subject to the equations y10 = y2 and y20 = −y1 ,together with initial values y1 (0) = 1, y2 (0) = 0. Then the function y1 (x) will be the cosinefunction.To use the routines above, we need to program the function f = f (x, y) that gives theright hand sides of the input ODE’s.
This is done as follows:> f:=proc(x,u,j)> if j=1 then RETURN(u[2]) else RETURN(-u[1]) fi:> end:That’s all. To run the programs we type the line> trapglobal(f,[1,0],.031415927,0,6.3,.0001,50):This means that we want to solve the system whose right hand sides are as given by f ,with initial condition array [1, 0]. The initial choice of the mesh size h is π/100, in order tofacilitate the comparison of the values that will be output with those of the cosine functionat the same points.
The range of x over which we are asking for the solution is from x = 0to x = 6.3. Our convergence criterion eps is set to .0001, and we are asking the programto print every 50th line of putput that it computes. Maple responds to this call with thefollowing output.At x=0h was reduced to.1570796350e-1.7853981750, [ .6845520546, -.7289635418]1.570796368, [-.0313962454, -.9995064533]2.356194568, [-.7289550591, -.6845604434]3.141592768, [-.9995063863, .0313843153]3.926990968, [-.6845688318, .7289465759]4.712389168, [ .0313723853,.9995063195]5.497787368, [ .7289380928,.6845772205]6.283185568, [ .9995062524, -.0313604551]68The Numerical Solution of Differential EquationsWe observe that the program first decided that the value of h that we gave it was too bigand it was cut in half.
Next we see that the accuracy is pretty good. At x = π we have 3or 4 correct digits, and we still have them at x = 2π. The values of the negative of the sinefunction, which are displayed above as the second column of unknowns, are less accurate atthose points. To improve the accuracy we might try reducing the error tolerance eps, butrealistically we will have to confess that the major source of imprecision lies in the Eulerand Trapezoidal combination itself, which, although it provides a good introduction to thephilosophy of these methods, is too crude to yield results of great accuracy over a long rangeof integration.2.11.2Example: The moon rocket in one dimensionAs a second example we will run the moon rocket in one dimension.
The equations thatwe’re solving now are given by (2.10.19). So all we need to do is to program the right handsides f , which we do as follows.> f:=proc(x,u,j)> if j=1 then RETURN(u[2]) else RETURN(-1/u[1]^2+.012/(60-u[1])^2) fi:> end:Then we invoke our routines by the statementtrapglobal(f,[1,1.4142],.02,0,250,.0001,1000):√This means that we are solving our system with initial data (1, 2), with an initial mesh sizeof h = .02, integrating over the range of time from t = 0 until t = 250, with a convergencetolerance eps of .0001, and printing the output every 1000 lines.
We inserted an extraline of program also, so as to halt the calculation as soon as y[1], the distance from earth,reached 59.75, which is the surface of the moon.Maple responded with the following output.At x=0h was reduced to.1000000000e-110.00000000, [7.911695068, .5027323920]20.00000000, [12.36222569, .4022135344]30.00000000, [16.11311118, .3523608113]40.00000000, [19.46812426, .3206445364]50.00000000, [22.55572225, .2979916199]60.00000000, [25.44558364, .2806820525]70.00000000, [28.18089771, .2668577906]80.00000000, [30.79081775, .2554698522]90.00000000, [33.29624630, .2458746369]100.0000000, [35.71287575, .2376534139]110.0000000, [38.05293965, .2305225677]120.0000000, [40.32630042, .2242857188]130.0000000, [42.54117736, .2188074327]140.0000000, [44.70467985, .2139997868]2.12 The big leaguesAt x=At x=At x=69150.0000000,160.0000000,170.0000000,174.2300000h was183.4500000h was188.0900000,211.7300000h was211.8100000,[46.82325670, .2098188816][48.90316649, .2062733168][50.95113244, .2034559376]increased to.2000000000e-1increased to.4000000000e-1[54.61091112, .2014073390]reduced to.2000000000e-1[59.75560109, .3622060968]So the trip was somewhat eventful.
The mesh size was reduced to .01 right away. It wasincreased again after 174 time units because at that time the rocket was moving quiteslowly, and again after 183 time units for the same reason. Impact on the surface of themoon occurred at time 211.81. Since each time unit corresponds to 13.5 minutes, this meansthat the whole trip took 211.8 · 13.5 minutes, or 47.65 hours – nearly two days.The reader might enjoy experimenting with this situation a little bit. For instance, ifwe reduce the initial velocity from 1.4142 by a small amount, then the rocket will reachsome maximum distance from Earth and will fall back to the ground without ever havingreached the moon.2.12The big leaguesThe three integration rules that we have studied are able to handle small-to-medium sizedproblems in reasonable time, and with good accuracy.
Some problems, however, are moredemanding than that. Our two-dimensional moon shot is an example of such a situation,where even a good method like the trapezoidal rule is unable to give the pinpoint accuracythat is needed. In this section we discuss a general family of methods, of which all threeof the rules that we have studied are examples, that includes methods of arbitrarily highaccuracy. These are the linear multistep methods.A general linear multistep method is of the formyn+1 =pXi=0a−i yn−i + hpX0b−i yn−i.(2.12.1)i=−1In order to compute the next value of y, namely yn+1 , we need to store p + 1 backwardsvalues of y and p + 1 backwards values of y 0 .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
