#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: