A look at discretization errors and irreversibility of numerical solvers as detailed in section 2.6 of the Borrelli and Coleman text.

Load the standard packages:

> map(with, {DEtools,plots}):

This is exercise 3: Solve y' + 2y = cos(t) and show all solutions converge to the particular solution 0.4cos(t) + 0.2sin(t), as t goes to infinity. dsolve would work fine, but this ODE is linear, so we'll proceed by hand (sort of). The integrating factor is exp(2*t), and, on the right hand side we are faced with integrating exp(2*t)*cos(t), which is messy enough to call on Maple for help.

> ans1:=int(exp(2*t)*cos(t),t) + c;

Divide through by exp(2*t), since after integrating, the left hand side is y*exp(2*t):

> ans2:=ans1/exp(2*t);

From that, we can read off the general solution:

> gensol:= 0.4*cos(t) + 0.2*sin(t) + c*exp(-2*t);

Since exp(-2*t) goes to zero, as t goes to infinity, all solutions will approach 0.4cos(t) + 0.2sin(t):

> limit(c*exp(-2*t),t=infinity);

Here, we verify that this particular solution, is actually a solution:

> yp:= 0.4*cos(t) + 0.2*sin(t):
leftside:=diff(yp,t)+2*yp;

That's also the right side, so yp satisfies the differential equation.

Part b asks for a plot showing several solutions on the rectangle t = [0, 20] and y = [-1.5, 1.5], illustrating that "all" solutions approach yp as t gets large. Here is a plot of yp:

> ypplot:=plot(yp,t=0..20,y=-1.5..1.5,
color=blue, thickness=2,axes=boxed):

ypplot;

Now, we'll use DEplot with several initial conditions. I'll build up the initial conditions using some typing and some uses of the sequence command.

> inits:=[ [y(0)=-1.5],
seq([y(4*k)=1.5], k=0..4), #these run across the top
seq([y(2+4*k)=-1.5],k=0..3)]; #these run across the bottom

Here is the DEplot:

> yplot:=DEplot(diff(y(t),t)=-2*y(t)+cos(t),y(t),
t=0..20,inits, y=-1.5..1.5,linecolor=red,
axes=boxed,arrows=none):

yplot;

That's pretty jagged, try a smaller stepsize and add some colors:

> yplot:=DEplot(diff(y(t),t)=-2*y(t)+cos(t),y(t),
t=0..22,inits, y=-1.5..1.5,
linecolor=[red,pink,black,green,yellow,violet,cyan,
aquamarine,plum,brown],
axes=boxed,arrows=none,stepsize=0.1):

yplot;

By the time t hits 19 or so, all those solutions are very close together. Here's the DEplot and the plot of yp together:

> display([ypplot, yplot],
view=[0..20,-1.5..1.5]);

One can imagine that a tiny error in estimating yp(20) would make any of the solutions shown above a possibility (if you were off by just a tiny bit, you might pick a point on any one of the different colored solutions). Picking a point on the brown solution, would give a very different value for y(10), than if you picked a point on the yellow solution! That is one thing which can go bad with numerical solvers. If many solutions, with radically different "beginnings", all come together near their "ends", starting at a point near the end and "backing up", can give bad results!

Here, we use rk4 with h = 0.1, and y(0) = 0.4 (the value of yp(0) ), and approximate y(20):

> rk4soln:=dsolve({diff(y(t),t)=-2*y(t)+cos(t), y(0)=0.4},
y(t), type=numeric,method=classical[rk4],
stepsize =0.1);

> rk4soln(20);

Now we solve backwards, using the point above, [20, y(20)] as an initial condition. If all goes well, when we evaluate this solution at t = 0, we should get y = 0.4.

> rk4back:=dsolve({diff(z(t),t)=-2*z(t)+cos(t), z(20)=.3458202008945779},
z(t), type=numeric,method=classical[rk4],
stepsize =0.1);

See how close z(0) is to 0.4?

> rk4back(0);

That's only off by 4 billion or so! I'll leave it to you to re-do the problem with h = 0.01 and h = 0.001.

>