> 
# 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, 29/Jan/96
# Math 308-509
# Chapter 8: Numerical Methods (for solving ODEs)
# Improved Euler's Method [ T-401/Eq (5) ]
# 
--------------------------------------------------------------------------------
# 
> 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];\
n:=convert((b-a)/h, rational);\
f:=(t, y)->1-t+4*y;\
deq:={diff(y(t), t) = f(t, y(t))};\
impeulerpts:=impeuler(f, i, h, n);\
list1:=convert(impeulerpts,listlist):\
plot(list1);

                                     a := 0

                                     b := 1

                                     h := .1

                                   i := [0, 1]

                                     n := 10

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

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

                   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]
                                  ])
--------------------------------------------------------------------------------
# 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];\
n:=convert((b-a)/h, rational);\
f:=(t, y)->1-t+4*y;\
deq:={diff(y(t), t) = f(t, y(t))};\
impeulerpts:=impeuler(f, i, h, n):\
for i from 0 to n by 2 do\
     impeulerpts[i]\
od;\
list2:=convert(impeulerpts,listlist):\
plot(list2);

                                     a := 0

                                     b := 1

                                    h := .05

                                   i := [0, 1]

                                     n := 20

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

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

                                     [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]
--------------------------------------------------------------------------------
# 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;\
deq:={diff(y(t), t) = f(t, y(t))};\
IC:={y(0)=1}; IVP:=deq union IC;\
sol:=dsolve(IVP, y(t));\
y:=unapply(subs(sol, y(t)), t);\
check:=IVP;\
for i from 0.0 to 1.0 by 0.1 do\
     [i, y(i)]\
od;\
plot(y(t), t=0..1);

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

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

                                IC := {y(0) = 1}

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

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

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

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

                                     [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({ list1,list2, y(t) }, t=0..1 );
--------------------------------------------------------------------------------
>