#Last update: Fri Aug 18 07:16:33 MDT 2006

#Wed Mar  6 11:59:18 2002
#Dr. Gustafson,
#I am in your 2250-4 class and am continually struggling with the codes for
#the problems in section 2.4, 2.5, and 2.6.  Actually 2.5 (Improved Euler)
#seems to be going pretty well, but I can't get very good  results for 2.4
#(Euler) or 2.6 (Runge-Kutta Idea).  Is there a website where I may be able
#to find help and/or codes on these sections, as nothing I type in for my
#codes will seem to work?  Thanks for your time.

# =========================================
# Can't copy with the mouse? A work-around:
# ========================================
# Run the application "xclipboard &" to capture the mouse copies of
# this file. Keep xclipboard near the xmaple window. Go to the netscape
# window, copy with the mouse. Switch to the xclipboard window. Copy
# with the mouse from the xclipboard. Then paste with mouse button 2 or
# mouse button 3 into xmaple.

# =========================================
# Maple code doesn't work? Read this:
# ========================================
# To type in a group, hold shift then press return, except
# for the last line of group, in which just return is used.
# If you copy multiple groups with the mouse, then split
# them using key F3 with the cursor placed at the front of
# a line where the split is to happen.

# Updated: Mon Sep 22 12:26:03 MDT 2003
# Warning: These snips of code made for y'=1-x-y, y(0)=3.
#          Code computes approx values for y(0.1) to y(1.0).
# 'Dots' is the list of dots for connect-the-dots graphics. 
# ========================================
# Euler. Group 1, initialize.
 f:=(x,y)->1-x-y:
 x0:=0:y0:=3:h:=0.1:Dots:=[x0,y0]:
# Group 2, repeat 10 times. Euler's method
 Y:=y0+h*f(x0,y0);
 x0:=x0+h:y0:=Y:Dots:=Dots,[x0,y0];
# Group 3, plot.
 plot([Dots]);

# ========================================
# Heun. Group 1, initialize.
 f:=(x,y)->1-x-y:
 x0:=0:y0:=3:h:=0.1:Dots:=[x0,y0]:
# Group 2, repeat 10 times. Heun method.
 Y1:=y0+h*f(x0,y0);
 Y:=y0+h*(f(x0,y0)+f(x0+h,Y1))/2:
 x0:=x0+h:y0:=Y:Dots:=Dots,[x0,y0];
# Group 3, plot.
 plot([Dots]);

# ========================================
# RK4. Group 1, initialize.
 f:=(x,y)->1-x-y:
 x0:=0:y0:=3:h:=0.1:Dots:=[x0,y0]:
# Group 2, repeat 10 times. RK4 method.
k1:=h*f(x0,y0):
k2:=h*f(x0+h/2,y0+k1/2):
k3:=h*f(x0+h/2,y0+k2/2):
k4:=h*f(x0+h,y0+k3):
Y:=y0+(k1+2*k2+2*k3+k4)/6:
x0:=x0+h:y0:=Y:Dots:=Dots,[x0,y0];
# Group 3, plot.
 plot([Dots]);

# Code snips for exact/error reports 2.4-#5
# =========================================
# Making multiple curves on one plot
# ========================================
  Exact:=x->2+x-exp(x);               # An exact solution for #5
  plot({Exact(x),[Dots]},x=0..1/2);   # plot exact and approx solutions

# ========================================
# How to create a Dots table for the exact solution
# ========================================
Exact:= x -> 2+x-exp(x):
ExactDots:=seq([Dots[j][1],Exact(Dots[j][1])],j=1..11);

# ========================================
# How to define and print percentage relative error:
# ========================================
  PRE:=unapply(evalf(100*abs(exact-approx)/abs(exact)),(exact,approx));
  ExactVal:=ExactDots[11][2]:   # Pick off exact y-value for x=0.5
  ApproxVal:=Dots[11][2]:       # Get Euler approx y-value for x=0.5
  PRE(ExactVal),ApproxVal);     # print percent relative error

# ========================================
# How to create a Dots table for percentage error
# ========================================
PRE:=unapply(evalf(100*abs(exact-approx)/abs(exact)),(exact,approx));
PREdots:=seq([Dots[j][1],PRE(ExactDots[j][2],Dots[j][2])],j=1..11);

# =========================================
# Printing results and tables
# Make tables with a pencil, it saves time. 
# ========================================
# To extract and print items 1,101,201,1001 from a list:
Dots1:=Dots[1],Dots[101],Dots[201],Dots[1001];

# =========================================
# Loop control
# ========================================
# To automate the production of the Dots list,
# add 1 below at the top of group 2, and 2 below
# to the end of group 2. 
#  1. for k from 1 to 10 do
#  2. od:
#
# Keyword "od:" is short for "end do:" -- use ":" to 
# stop loop results from printing.
# the end