> 
# This worksheet demonstrates the use of the "Share" library routine 
# rungekuttahf, a hardware floating point version of the classical
# Runge Kutta method. It graphs the Runge Kutta method solution 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)
# Runge-Kutta Method (classical) [ T-407/Eqs 2-3 ]
# -> hardware floating point chip used
#      (VERY FAST: "It rules.")
# 
> restart;
> with(share): readshare(ODE,plots):
See ?share and ?share,contents for information about the share library
> ?rungekuttahf
# 
# T-408/Table 8.4.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.2; i:=[0, 1];

                                a := 0


                                b := 1


                               h := .2


                             i := [0, 1]

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

                                n := 5

> 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

> rkhfpts:=rungekuttahf(f, i, h, n);

rkhfpts := array(0 .. 5, [

    (0) = [0, 1.]

    (1) = [.2000000000000000, 2.501600000000000]

    (2) = [.4000000000000000, 5.777635840000000]

    (3) = [.6000000000000001, 12.99717789081600]

    (4) = [.8000000000000000, 28.98076814454948]

    (5) = [1., 64.44157912444676]

    ])

> plot(makelist(rkhfpts));

# T-403/Table 8.3.1: Column 4, via raw dump; graph of solution.
> 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

> rkhfpts:=rungekuttahf(f, i, h, n);

rkhfpts := array(0 .. 10, [

    (0) = [0, 1.]

    (1) = [.1000000000000000, 1.608933333333333]

    (2) = [.2000000000000000, 2.505006151111111]

    (3) = [.3000000000000000, 3.829414509150815]

    (4) = [.4000000000000000, 5.792785270450575]

    (5) = [.5000000000000000, 8.709317547440138]

    (6) = [.6000000000000000, 13.04771262943470]

    (7) = [.7000000000000000, 19.50714785308206]

    (8) = [.7999999999999999, 29.13060935737094]

    (9) = [.8999999999999999, 43.47395433203548]

    (10) = [1.000000000000000, 64.85810680890840]

    ])

> plot(makelist(rkhfpts));

# ASIDE: Every other item in rkpts:
> for i from 0 to n by 2 do
>      rkhfpts[i]
> od;

                               [0, 1.]


                [.2000000000000000, 2.505006151111111]


                [.4000000000000000, 5.792785270450575]


                [.6000000000000000, 13.04771262943470]


                [.7999999999999999, 29.13060935737094]


                [1.000000000000000, 64.85810680890840]

# T-408/Table 8.4.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 = 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);

# END OF TRANSMISSION.
> #STOP!
# 
>