# This worksheet demonstrates the use of the "Share" library routine 
# firsteuler.  It graphs the 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)
# Euler's Method [ T-389/Eq (6) ]
# 
# 
> with(share): readshare(ODE,plots):
See ?share and ?share,contents for information about the share library
> ?firsteuler
# 
# T-390/Table 8.1.1: Column 1, 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

> eulerpts:=firsteuler(f, i, h, n);

eulerpts := array(0 .. 10, [

    (0) = [0, 1.]

    (1) = [.1, 1.5]

    (2) = [.2, 2.19]

    (3) = [.3, 3.146]

    (4) = [.4, 4.4744]

    (5) = [.5, 6.32416]

    (6) = [.6, 8.903824]

    (7) = [.7, 12.5053536]

    (8) = [.8, 17.53749504]

    (9) = [.9, 24.57249306]

    (10) = [1.0, 34.41149028]

    ])

> plot(makelist(eulerpts));

# T-390/Table 8.1.1: Column 2, 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

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

                               [0, 1.]


                            [.10, 1.5475]


                          [.20, 2.32490000]


                          [.30, 3.433356000]


                          [.40, 5.018532640]


                          [.50, 7.290187002]


                          [.60, 10.55036928]


                          [.70, 15.23403177]


                          [.80, 21.96750574]


                          [.90, 31.65270827]


                         [1.00, 45.58839992]

> plot(makelist(eulerpts));

# T-390/Table 8.1.1: Column 3, via judicious paring (every 4th item);
# graph of solution.
> unassign('t', 'y');
> a:=0; b:=1; h:=0.025; i:=[0, 1];

                                a := 0


                                b := 1


                              h := .025


                             i := [0, 1]

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

                               n := 40

> 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

> eulerpts:=firsteuler(f, i, h, n):
> for i from 0 to n by 4 do
>      eulerpts[i]
> od;

                               [0, 1.]


                         [.100, 1.576118750]


                         [.200, 2.408011713]


                         [.300, 3.614383699]


                         [.400, 5.369030425]


                         [.500, 7.926406197]


                         [.600, 11.65905757]


                         [.700, 17.11242994]


                         [.800, 25.08510991]


                         [.900, 36.74630817]


                         [1.000, 53.80786605]

> plot(makelist(eulerpts));

# T-390/Table 8.1.1: Column 4, via judicious paring (every 10th item);
# graph of solution.
> unassign('t', 'y');
> a:=0; b:=1; h:=0.01; i:=[0, 1];

                                a := 0


                                b := 1


                               h := .01


                             i := [0, 1]

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

                               n := 100

> 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

> eulerpts:=firsteuler(f, i, h, n):
> for i from 0 to n by 10 do
>      eulerpts[i]
> od;

                               [0, 1.]


                          [.10, 1.595290089]


                          [.20, 2.464458735]


                          [.30, 3.739034547]


                          [.40, 5.613712001]


                          [.50, 8.376686480]


                          [.60, 12.45455756]


                          [.70, 18.47879679]


                          [.80, 27.38413639]


                          [.90, 40.55420836]


                         [1.00, 60.03712599]

> plot(makelist(eulerpts));

# T-390/Table 8.1.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(0) = 1, -- y(t) = 1 - t + 4 y(t)}
                               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/4 + 19/4 exp(4 t) = 1/4 + 19/4 exp(4 t), 1 = 1}

> 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);

# END OF TRANSMISSION.
> #STOP!
>