> 
# This worksheet demonstrates the use of the "Share" library routine 
# Runge Kutta.  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, 29/Jan/96
# Math 308-509
# Chapter 8: Numerical Methods (for solving ODEs)
# Runge-Kutta Method (classical) [ T-407/Eqs 2-3 ]
# 
--------------------------------------------------------------------------------
# 
> with(share):\
readshare(ODE,plots):
See ?share and ?share,contents for information about the share library

> ?rungekutta
--------------------------------------------------------------------------------
# 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];\
n:=convert((b-a)/h, rational);\
f:=(t, y)->1-t+4*y;\
deq:={diff(y(t), t) = f(t, y(t))};\
rkpts:=rungekutta(f, i, h, n);\
plot(makelist(rkpts));
--------------------------------------------------------------------------------
# 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];\
n:=convert((b-a)/h, rational);\
f:=(t, y)->1-t+4*y;\
deq:={diff(y(t), t) = f(t, y(t))};\
rkpts:=rungekutta(f, i, h, n);\
plot(makelist(rkpts));
--------------------------------------------------------------------------------
# ASIDE: Every other item in rkpts:
> for i from 0 to n by 2 do\
     rkpts[i]\
od;
--------------------------------------------------------------------------------
# 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;\
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);
--------------------------------------------------------------------------------
# END OF TRANSMISSION.
> #STOP!
--------------------------------------------------------------------------------
# 
>