# # The idea for higher order versions of Euler is introduced. Maple # procedures (3) are provided to construct this new method. The solution # approximated by this technique is compared to the exact solution. # The trade-offs are mentioned for this type of method. # # Written by Doug Hensley # # HIGHER ORDER TAYLOR'S SERIES METHODS # The idea here is that if we have a DE such as y'=xy, we can differentiate it. # Taking the derivative of this DE gives another DE, this time involving y'': # y''=(xy)'=y+x(y'), and this=y+x(xy) because y'=xy. Now if the basic Euler # method makes steps based on the linear extrapolation dy==slope*dx, # slope=y' =x_ky_k, then a second-order Taylor's version of the Euler method # would extrapolate via dy==slope*dx+second*dx^2/2!, where # "second"=y+xy@x=x_k,y=y_k. The code and output below give an # implementation of this idea, with implementation carried out on the # (separable)DE y'=xy. The code is illustrated on a problem for which it is not # needed so that the answer the code gives can be compared with the real # answer! > fn:=(x,y)->x*y; fn := (x,y) -> x y -------------------------------------------------------------------------------- # # The arcana of this code is the use of the device deriv.j which constructs a # sequence of distinct names for the functions which are to be stored in the list # of derivatives. This CAS is balky when it comes to doing this sort of work, # and some natural-seeming approaches run aground on the hidden reefs of # Maple. > derivmaker:=proc(f,depth)local derivlist,deriv,j,k,s,t;\ deriv.0:=f;derivlist:=f;for k from 1 to depth do j:=k-1;deriv.k:=unapply(diff(deriv.j(s,t),s)+f(s,t)*diff(deriv.j(s,t),t),(s,t)) od;[seq(deriv.k,k=0..depth)]end;\ derivmaker := proc(f,depth) local derivlist,deriv,j,k,s,t; deriv.0 := f; derivlist := f; for k to depth do j := k-1; deriv.k := unapply(diff(deriv.j(s,t),s)+f(s,t)*diff(deriv.j(s,t),t),s,t) od; [seq(deriv.k,k = 0 .. depth)] end # This code takes as input a function of two variables, and a depth...how many # derivatives do you want?. In hand calculation at the top, you saw depth=2. # )The code returns a list of functions of two variables. -------------------------------------------------------------------------------- > applyem:=proc(f,x,y)f(x,y)end; applyem := proc(f,x,y) f(x,y) end > # # This code takes as input a function of two variables, and two numbers or # expressions, and returns the expression got by setting the two variables equal # to those two expressions. This is a fine distinction. The object (x,y)->x*y is a # function. It is the same function as the object (s,t)->s*t, and neither of these # actually have anything to do with x,y,s or t. The expression x*y, however, is # not a function (not to Maple, anyhow, and s*t is yet another expression. The # point of keeping the derivatives as functions in the first place is that this way, # we're not tied to any specific choice of variable names. We can later use this # code on, say, u'=f(u,t) and it will work. -------------------------------------------------------------------------------- > newymaker:=proc(taylist,x,y,h)local tot,k,n,hk;n:=nops(taylist);tot:=y;hk:=h;\ for k from 1 to n do tot:=evalf(tot+hk*applyem(taylist[k],x,y)/k!);hk:=hk*h od;tot end;\ newymaker := proc(taylist,x,y,h) local tot,k,n,hk; n := nops(taylist); tot := y; hk := h; for k to n do tot := evalf(tot+hk*applyem(taylist[k],x,y)/k!); hk := hk*h od; tot end > # # This code takes as input the list of derivatives (as functions), the current # values of x and y, and the step size, and returns the Taylor's polynomial # approximation y+a_1h+a_2h^2+..., where a_1=y'@(x,y), # a_2=(1/2!)y''@(x,y) etc. and the calculation is carried out to the number of # terms specified by "depth". -------------------------------------------------------------------------------- > # This code takes as input the function (as, for instance, (x,y)->x*y ), the # starting and ending values of x, the starting value of y, the step size, and the # depth. It returns the value of y corresponding to the ending value of x, and as # a global spinoff, the list (ptlist) of the points calculated along the way. > MyTaylorDE:=proc(f,firstx,lastx,firsty,stepsize,depth)local x,y,s,t,u,v;x:=firstx;y:=firsty;s:=NULL;\ t:=derivmaker(f,depth);\ while x MyTaylorDE(fn,1,4,1,1,3); [[2, 4.083333334], [3, 39.30208335], [4, 815.5182299]] -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- > MyTaylorDE((x,y)->x*y,1,4,1,1,3); [[2, 4.083333334], [3, 39.30208335], [4, 815.5182299]] -------------------------------------------------------------------------------- > MyTaylorDE((x,y)->x*y,1,4,1,1/2,3); [[3/2, 1.859375000], [2, 4.418739319], [5/2, 13.38279643], [3, 51.45972749], [7/2, 250.0621133], [4, 1527.438322]] -------------------------------------------------------------------------------- > MyTaylorDE((x,y)->x*y,1,4,1,1/4,3); [[5/4, 1.324544271], [3/2, 1.867397532], [7/4, 2.802179079], [2, 4.475324407], [9/4, 7.606740363], [5/2, 13.75905551], [11/4, 26.48251480], [3, 54.23362313], [13/4, 118.1593879], [7/2, 273.8434474], [15/4, 675.0063811], [4, 1769.351567]] -------------------------------------------------------------------------------- -------------------------------------------------------------------------------- > MyTaylorDE((x,y)->x*y,1,4,1,1/8,3); -------------------------------------------------------------------------------- > MyTaylorDE((x,y)->x*y,1,4,1,1/8,4); -------------------------------------------------------------------------------- > MyTaylorDE((x,y)->x*y,1,4,1,1/8,5); -------------------------------------------------------------------------------- > MyTaylorDE((x,y)->x*y,1,4,1,1/8,6); -------------------------------------------------------------------------------- > MyTaylorDE((x,y)->x*y,1,4,1,1/1000,3); -------------------------------------------------------------------------------- > MyTaylorDE((x,y)->x*y,1,4,1,1/8,7); -------------------------------------------------------------------------------- # The exact solution of this DE is, in general, y=Ae^(x^2/2). For the specific # initial conditions x=1,y=1, the solution is y=e^((x^2-1)/2). > evalf(E^(15/2)); -------------------------------------------------------------------------------- # Our two last approximations bracket the exact solution. The first of these had # a great many steps, but a relatively low-order approximation. The second # involved a high-order extrapolation, but a seven-league boots step size of # 1/8. Note the trade-off: the latter method involves more abstract calculation, # with the list below representing part of the abstract overhead, before overt # arithmetic ever gets going. On the other hand, when it finally does come to # arithmetic, there is a lot more of it when working with a lower-order # approximation! The better computers get, the stronger this effect will get. > > listderivs:=derivmaker(fn,4):listinxy:=map(applyem,listderivs,x,y):map(expand,listinxy); -------------------------------------------------------------------------------- --------------------------------------------------------------------------------