#Last update: December 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? Work-arounds:
# ========================================
# Try first: hold down SHIFT and then hit the paste button on
# the mouse. Usually needed from xterm windows and firefox to other
# edit windows and browser forms.
#
# If pasting into Xmaple, then try just the mouse paste button. It
# might work. For some systems, ctrl-V will perform the paste.
#
# Lastly: Run the application "xclipboard &" to capture mouse copies.
# Keep xclipboard near the xmaple window. Go to the browser 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 case use just return.
# 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: December 2006
# 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).
# '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]:n:=10:
# Group 2, repeat n times. Euler's method
for i from 1 to n do
Y:=y0+h*f(x0,y0);
x0:=x0+h:y0:=Y:Dots:=Dots,[x0,y0];
od:
# Group 3, display dots and plot.
Dots;
plot([Dots]);
# ========================================
# Heun. Group 1, initialize.
f:=(x,y)->1-x-y:
x0:=0:y0:=3:h:=0.1:Dots:=[x0,y0]:n:=10:
# Group 2, repeat n times. Heun method.
for i from 1 to n do
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];
od:
# Group 3, display dots and plot.
Dots;
plot([Dots]);
# ========================================
# RK4. Group 1, initialize.
f:=(x,y)->1-x-y:
x0:=0:y0:=3:h:=0.1:Dots:=[x0,y0]:n:=10:
# Group 2, repeat n times. RK4 method.
for i from 1 to n do
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];
od:
# Group 3, display some dots and plot.
Dots[1],Dots[2],Dots[n+1];
plot([Dots]);
# Code snips for exact/error reports
# =========================================
# Making multiple curves on one plot
# ========================================
Exact:=x->2-x+exp(-x); # An exact solution
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):n:=10:
ExactDots:=seq([Dots[j][1],Exact(Dots[j][1])],j=1..n+1);
# ========================================
# How to define and print percentage relative error:
# ========================================
P:=unapply(evalf(100*abs(exact-approx)/abs(exact)),(exact,approx));
ExactVal:=Exact(Dots[11][1]): # Compute exact y-value for x=1
ApproxVal:=Dots[11][2]: # Get Euler approx y-value for x=1
P(ExactVal),ApproxVal); # print percent relative error
# ========================================
# How to create a Dots table for percentage error
# ========================================
P:=unapply(evalf(100*abs(exact-approx)/abs(exact)),(exact,approx));
Pdots:=seq([Dots[j][1],P(Exact(Dots[j][1]),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 a Dots list,
# enclose the desired code between 1 and 2 below.
# 1. for k from 1 to 10 do
# 2. od:
#
# Keyword "od:" is short for "end do:"
# Use ":" to stop loop results from printing.
# =========================================
# Debug
# ========================================
# To remove loop control and do it by hand, insert
# pound (#) signs as follows:
# 1. # for k from 1 to 10 do
# 2. # od:
# The hand-done loop is made by placing the mouse cursor
# in the group, then press return. Repeat for each loop step,
# which is 10 times for the loop above.
# the end