#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