> 
# This worksheet demonstrates the use of the "Share" library routine 
# impeuler.  It graphs the Improved Euler method solution for h = 
# .1,.05,.025, and .01 along with the actual solution from dsolve.  It 
# then plots all of these together to show a comparision. 
#
# Art Belmonte
# Mon, 27/May/96
# Math 308-509 [Maple V Release 4]
# 
# Chapter 8: Numerical Methods (for solving ODEs)
# Improved Euler's Method [ T-401/Eq (5) ]
# 
# Math 308 autoloading (*ALL* current/future routines)
# (DO ONCE PER SEMESTER):
# cp  ~belmonte/.mapleinit  ~
# 
# On-line help; files also in
# /courses/math308/INSTRUCTORS/belmonte/DOCS
# 
> restart;
> with(share): readshare(ODE,plots):
See ?share and ?share,contents for information about the share library
> ?impeuler
# 
# T-403/Table 8.3.1: Column 3, via raw dump; graph of solution (all
# plots omitted to save space).
# LEGEND:
#      a - left endpoint of interval for independent variable
#      b - right endpoint of interval for independent variable
#      f  - derivative function in D.E.; from dy/dt = f(t, y)
#      h - step size
#      i  - list of initial values: [ t0, y0 ] = [ indep_var, dep_var ]
#      n - number of steps: (b-a)/h, as a RATIONAL!
> unassign('t', 'y');
> a:=0; b:=1; h:=0.1; i:=[0, 1];

                                a := 0


                                b := 1


                               h := .1


                             i := [0, 1]

> n:=convert((b-a)/h, rational);

                               n := 10

> f:=(t, y)->1-t+4*y;

                      f := (t, y) -> 1 - t + 4 y

> deq:={diff(y(t), t) = f(t, y(t))};

                          d
                  deq := {-- y(t) = 1 - t + 4 y(t)}
                          dt

> impeulerpts:=impeuler(f, i, h, n);

impeulerpts := array(0 .. 10, [

    (0) = [0, 1.]

    (1) = [.1, 1.595000000]

    (2) = [.2, 2.463600000]

    (3) = [.3, 3.737128000]

    (4) = [.4, 5.609949440]

    (5) = [.5, 8.369725171]

    (6) = [.6, 12.44219325]

    (7) = [.7, 18.45744601]

    (8) = [.8, 27.34802009]

    (9) = [.9, 40.49406974]

    (10) = [1.0, 59.93822322]

    ])

> plot(makelist(impeulerpts));

# T-403/Table 8.3.1: Column 4, via judicious paring (every 2nd item);
# graph of solution.
> unassign('t', 'y');
> a:=0; b:=1; h:=0.05; i:=[0, 1];

                                a := 0


                                b := 1


                               h := .05


                             i := [0, 1]

> n:=convert((b-a)/h, rational);

                               n := 20

> f:=(t, y)->1-t+4*y;

                      f := (t, y) -> 1 - t + 4 y

> deq:={diff(y(t), t) = f(t, y(t))};

                          d
                  deq := {-- y(t) = 1 - t + 4 y(t)}
                          dt

> impeulerpts:=impeuler(f, i, h, n):
> for i from 0 to n by 2 do
>      impeulerpts[i]
> od;

                               [0, 1.]


                          [.10, 1.604975000]


                          [.20, 2.493209790]


                          [.30, 3.803048452]


                          [.40, 5.740402317]


                          [.50, 8.611749809]


                          [.60, 12.87325342]


                          [.70, 19.20386539]


                          [.80, 28.61413826]


                          [.90, 42.60817839]


                         [1.00, 63.42469773]

> plot(makelist(impeulerpts));

# T-403/Table 8.3.1: Column 5, exact solution (via dsolve); graph of
# solution.
> unassign('t', 'y');
> f:=(t, y)->1-t+4*y;

                      f := (t, y) -> 1 - t + 4 y

> deq:={diff(y(t), t) = f(t, y(t))};

                          d
                  deq := {-- y(t) = 1 - t + 4 y(t)}
                          dt

> IC:={y(0)=1}; IVP:=deq union IC;

                           IC := {y(0) = 1}


                     d
             IVP := {-- y(t) = 1 - t + 4 y(t), y(0) = 1}
                     dt

> sol:=dsolve(IVP, y(t));

                                             19
              sol := y(t) = - 3/16 + 1/4 t + -- exp(4 t)
                                             16

> y:=unapply(subs(sol, y(t)), t);

                                           19
                y := t -> - 3/16 + 1/4 t + -- exp(4 t)
                                           16

> check:=IVP;

     check := {1 = 1, 1/4 + 19/4 exp(4 t) = 1/4 + 19/4 exp(4 t)}

> for i from 0.0 to 1.0 by 0.1 do
>      [i, y(i)]
> od;

                                [0, 1]


                          [.1, 1.609041829]


                          [.2, 2.505329852]


                          [.3, 3.830138846]


                          [.4, 5.794226004]


                          [.5, 8.712004118]


                          [.6, 13.05252195]


                          [.7, 19.51551804]


                          [.8, 29.14487961]


                          [.9, 43.49790340]


                          [1.0, 64.89780316]

> plot(y(t), t=0..1);