#Header Begins--Do not remove##################################################
#
# This file suggests the format of some routines that will help
# you to write forward differencing code.
#
#
# (c) Eugene Fiume, University of Toronto, 1995.
#
# You are permitted to use this code for noncommercial, educational
# purposes only. This header must not be removed from this file.
#
#Header Ends--Do not remove####################################################
with(plots):
#
# Modify this value of epsilon accordingly to investigate the accuracy of
# forward differencing.
#
eps := 0.0000000001:
eps := 0.0:
#
# Define a polynomial in indeterminate t given coefficients A as in:
# [t^(n-1) t^(n-2) ... t 1] dot A, where n=nops(A).
#
poly := proc(A,t)
local i,n;
n := nops(A);
sum(A[i]*t^(n-i), i=1..n);
end:
#
# Compute forward differences for monomial
# [t^(n-1) t^(n-2) ... t 1] dot A, where n=nops(A).
#
# Return in the form:
# [p(t), Dp(t), DDp(t), DDDp(t), ... ]
#
# where "DDD" means the nth forward difference (in this case, 3).
#
# Sample output:
#
# > fwddiff([a,b,c,d],t,h);
# 3 2 2 2 2 3
# [a t + b t + c t + d, 3 a h t + (2 b h + 3 a h ) t + b h + c h + a h ,
#
# 2 3 2 3
# 6 a h t + 6 a h + 2 b h , 6 a h ]
#
# This should not be a surprising result.
#
fwddiff := proc (A,t,dt)
end:
#
# Initialise the forward differences of polynomials P (created by fwddiff())
# in an indeterminate t to t0, and with dt set accordingly to specific value
# dt0. Note: normally h is the indeterminate specifying the step size, and
# dt0 is the specific value for the indeterminate, i.e., one divided by the
# number of samples (less one) that you want to take.
#
# Example:
# > initdiffs(fwddiff([a,b,c,d],t,h), t, 0, h, h); #keeps the result symbolic
#
# 2 3 3 2 3
# [d, b h + c h + a h , 6 a h + 2 b h , 6 a h ]
#
# > initdiffs(fwddiff([a,b,c,d],t,h), t, 0, h, 1/16); # gives a value for h
#
# [d, 1/256 b + 1/16 c + 1/4096 a, 3/2048 a + 1/128 b, 3/2048 a]
#
#
initdiffs := proc(P, t, t0, dt, dt0)
end:
#
# Evaluate polynomial on k+1 equally spaced points on [t0,t1]
# using forward differencing. Produce points in the format
# [ [t0,p(t0)], [t1,p(t1)], ... , [tk, p(tk)]].
#
# You will need to give the polynomial's coefficients as a vector A
# (see above). The values t0 and t1 give the lower and upper bounds
# of the interval, and k gives the number of steps as above. You
# need to use initdiffs and fwddiff to create and initialise the
# forward differences. Look above for a *big* hint. Then you
# need to accumulate points on the curve iteratively throught the
# forward-differencing process (i.e., addition).
#
# Later, you will need to modify this routine to use eps above to
# simulate an error in the nth forward difference.
#
evalpoly := proc(A,t0,t1,k)
end:
#
# Invoke this routine to plot "exactly" computed polynomials against
# the forward differencing solution. A, t0, t1 and k are as above.
#
# Sample execution:
#
# plotpoly([1,2,3,4], 0, 1, 10);
#
# should give a gently upward sloping polynomial on [0,1].
#
plotpoly := proc(A,t0,t1,k)
local t,plot1,plot2,plot3, E;
E := evalpoly(A,t0,t1,k);
plot1 := plot(E, style=point, colour=GREEN);
plot2 := plot(E, style=line, colour=RED);
plot3 := plot(poly(A,t), t=t0..t1, colour=YELLOW);
display({plot1,plot2,plot3});
end: