Name:

It is critical that one executes the commands in this file, from top to bottom, as I have broken a Maple rule of thumb, which suggests that you never re-use names. Executing the commands out of order, will cause unpredictable results.

Load some packages:

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

Warning, new definition for translate

For Euler's Method, it is best to think of t - 3 (the "distance" from t to 3) as a quantity h. Moving h units to the right, gives a y value of 4 + (-2)*h. Euler's method is simply the process of moving h units, recording the new y value, using the new t value (3 + h) and the new y value to calculate the new slope, y'(3 + h) (using the differential equation), and then again moving h units to the right, etc. Take y' = t^2 / y, y(3) = 4. Since the forcing function f(t,y) = t^2 / y and its partial with respect to y, -t^2 / y^2 are continuous away from y = 0 (the x-axis), we have existence and uniqueness of solutions, if we avoid the x-axis. Using h = 0.2, let's estimate y(4), by hand. We'll build a list of pairs, [t, y], approximating points on the graph of y(t). Here is our start of that list of lists, the value of h, a function for the right hand side, f(t,y), and the "real solution":

> pairs:=[[3.,4.]];
h:=0.2;
f:=(t,y)->t^2/y;
ans:=dsolve({diff(y(t),t)=f(t,y(t)),y(3)=4},y(t));

Euler's method sets the next y value to be the y-value on the tangent line, h units to the right of the last t-value. That is:

> nexty:= 4+h*f(3.,4.);

and the next t value to be , where and are the last t and y-values respectively.

> nextt:=3.+h;

Now, add this new pair to the list of pairs of points (op(pairs), just removes the outermost list brackets of pairs, so this command says "make a list named pairs which is the old list, followed by the pair [nextt, nexty]"):

> pairs:=[ op(pairs),[nextt,nexty] ];

A (fancy) plot of what we've done. The true solution is red, the Euler numerical solution is blue, and the pairs are black boxes:

> plot([pairs,pairs,rhs(ans)],t=2.8..3.4,
style=[point,line,line],symbol=BOX,
color=[black,blue,red]);

The next pair by direct typing:

> nexty:=4.45+h*f(3.2, 4.45);
nextt:=3.2+h;

Add this new pair to the list of pairs:

> pairs:=[ op(pairs),[nextt,nexty] ];

And a plot:

> plot([pairs,pairs,rhs(ans)],t=2.8..3.6,
style=[point,line,line],symbol=BOX,
color=[black,blue,red]);

Now we simply "age" the pair [nextt, nexty], since they are now the old t and old y:

> oldt:=nextt;
oldy:=nexty;

and use Euler's method to get the next values:

> nexty:=oldy+h*f(oldt, oldy);
nextt:=oldt+h;
pairs:=[ op(pairs), [nextt, nexty] ];

> plot([pairs,pairs,rhs(ans)],t=2.8..3.8,
style=[point,line,line],symbol=BOX,
color=[black,blue,red]);

That was easy, this is faster:

> oldt:=nextt:
oldy:=nexty:
nexty:=oldy+h*f(oldt, oldy);
nextt:=oldt+h;
pairs:=[ op(pairs), [nextt, nexty] ];

> plot([pairs,pairs,rhs(ans)],t=2.8..4,
style=[point,line,line],symbol=BOX,
color=[black,blue,red]);

Perform the last iteration of Euler's method to estimate y(4). What is your estimate of y(4)? Show a plot of the real solution, and the 5 steps of your estimate.

>

Example 2.5.2, if y' = y*sin(3*t), and y(0) = 1, estimate y(4), can be handled quickly with the following loop, or repeated iteration of Euler's method outlined above. We start by clearing (unassigning) all important values:

> pairs:='pairs':
nextt:='nextt':
nexty:='nexty':
h:='h':
f:='f':
t:='t':

Next, we define the rhs function f, the list of pairs, i.e. just the empty list, the value of the stepsize h, and the first value of oldy:

> f:= (t,y) -> y*sin(3*t):
pairs:=[ [0,1] ]:
h:=0.1:
oldy:=1:

We use a Maple looping procedure "for", to run the steps introduced above, as t runs from 0 to 4, with steps of size h. Since t is automatically increased by h, each iteration of the loop, we can do away with the nextt calculation.

> for tval from 0 to 4-h by h do
nexty:=oldy+h*f(tval, oldy):
pairs:= [ op(pairs), [tval+h, nexty] ]:
oldy:=nexty:
od:

Now, that was a breeze!

Here, we make a matrix of the points in pairs:

> evalm([ [`t`, `y(t)`],op(pairs[1..13])]),
evalm( [ [`t`, `y(t)`],op(pairs[14..27])]),
evalm([ [`t`, `y(t)`],op(pairs[28..41])]);

>

Here's our fancy plot, after we solve the differential equation:

> ans:=dsolve({diff(y(t),t)=f(t,y(t)),y(0)=1},y(t)):

plot([pairs,pairs,rhs(ans)],t=-0.3..4.3,
style=[point,line,line],symbol=BOX,
color=[black,blue,red]);

Make a plot like the one above with h = 0.05

>

Here is a fast way to get Euler's method with a stepsize of 0.1. The option foreuler means forward Euler.

> deq1:=diff(y(t),t)=y(t)*sin(3*t):

init1:=y(0)=1:

ans1 := dsolve({deq1, init1}, y(t), type=numeric,
method=classical[foreuler], stepsize=0.1);

Evaluating this gives a list of equations:

> ans1(0);

To get pairs [t,y], we can use the substitute command:

> subs(ans1(0.1), [t, y(t)]);

We can also turn ans1 into a function, soln, which on input u gives y(u) via the Euler method:

> soln := u -> subs(ans1(u),y(t));

Now, we just use soln(0.2) to get a specific value, or soln as a plot argument

> soln(0.2);

> plot(soln,0..4,numpoints=41,adaptive=false);

We can also make a list of points and plot them as data:

> pairs:=[seq([h*k, soln(h*k)], k=0..40)];

> plot(pairs,style=point,symbol=BOX,color=blue);

In the command above, which defines ans1, if we replace method=classical[foreuler], with method=classical[heunform], Maple will use the Heun method, and method=classical[rk4] gives the Runge-Kutta method of order 4.

Do problems 1 and 3 of section 2.5

>