Math 640: EXPERIMENTAL MATHEMATICS (Symbolics Meets Numerics!)
Spring 2014 (Rutgers University) Webpage
http://sites.math.rutgers.edu/~zeilberg/math640_14.html
Last Update: May 7, 2014
 Teacher:
Dr. Doron ZEILBERGER ("Dr. Z")

Classroom:
Allison Road Classroom Building
[Busch Campus], IML Room 118

Time: Mondays and Thursdays , period 3 (12:00noon1:20pm)

"Textbooks": Links to wikipedia and other internet sites

Dr. Zeilberger's Office: Hill Center 704

Dr. Zeilberger's Email:
zeilberg at math dot rutgers dot edu (you MUST have MathIsFun in the message,
or ExpMathRocks)

Dr. Zeilberger's Office Hours (Spring 2014): MTh 10:30am11:30am and by appointment
Description
Experimental Mathematics used to be considered an oxymoron,
but the future of mathematics is in that direction.
In addition to learning the
philosophy and methodology of this budding field, students will
become computeralgebra wizards, and that should be very helpful in whatever
mathematical specialty they are doing (or will do) research in.
We will first learn Maple, and how to program in it.
This semester we will do Numerical Analysis (Numerical solutions of ODEs (RungeKutta)
and PDEs(using both finite difference schemes and the finite element method)) via Symbolic Computation. We will also experiment with
Approximation Theory, and see how symbolic computation can help it.
In addition to the actual, very important content, students will
master the methodology of computergenerated and
computerassisted research that is so crucial for their future.
There are no prerequisites, and no previous programming knowledge is assumed.
Also, no overlap with previous years.
The final projects will with high probability lead to publishable articles,
and with strictly positive probability, to seminal ones.
"Textbooks"
Added March 23, 2014:
Pick a project
Getting Started Homework (assigned Jan. 12, 2014):
Read
L. Nick Trefethen's Definition of Numerical Analysis.
and
L. Nick Trefethen's beautiful essay on Numerical Analysis, pp. 604615 from this url.
How to submit your homework
If you don't already have a public_html directory, go
to commandline Unix (or Linux), and type

mkdir public_html

chmod 755 public_html

cd public_html

mkdir em14

chmod 755 em14

cd em14

create the homework file, call it hw0.txt (and in the future, hw1.txt , hw2.txt etc.)
and do

chmod 755 *
More Getting Started Homework (assigned Jan. 21, 2014):
Create a file hw0.txt in the above em14 directory, call it hw0.txt
that says
"I love Experimental Mathematics"
Diary and Homework
Programs done on Thurs., Jan. 23, 2014
C1.txt ,
Contains procedures

M(a,b): A "homemade" program for multiplying two positive integers

K(a,b): A faster version using Karatsuba's
algorithm.
Homework for Thurs., Jan. 23, class (due Sunday Jan. 26, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw1.txt

(For novices only)
Read and do all the examples, plus make up similar ones,
in the first 30 pages of Frank Garvan's
awesome Maple booklet
(
part 1,
part 2)
Only record a small sample in hw1.txt .
Note: a few commands are no longer valid in today's Maple. The
most important one is that " has been replaced by %.

(For novices only) Extend procedures M(a,b) and K(a,b) of
C1.txt ,
to write procedures Mbase(a,b,B), and Kbase(a,b,B), that does things in base B.
Of course, since Maple uses base 10, the input and output are still the same,
but experiment with different B, and see whether you get the same answers.

(For everyone) Get ready for the next class by reading the wikipedia article on
Newton's method

(Mandatory for experts, optional for novices)
Implement Newton's method. Write a program NR(f,x,x0,eps) that inputs an expression f in the
variable x, a variable x, an initial guess x0, and a small positive number eps, and outputs
the root of f(x)=0 to error less than eps. For example
NR(x^22,x,1,10^(10));
should be the same as evalf(sqrt(2));

(Mandatory for experts, optional for novices)
Using the time(..) command,
write programs CompM(N), CompK(N), that inputs a positive integer N, and
outputs a list of length N, whose ith entry is the average running time
of M(a,b) (and K(a,b)) respectively
of ten randomlydrawn pairs of integers with d digits.
See whether they are asymptotically a constant times d^{2} and d^{(log[2](3))}=
d^{1.584962501}

(Challenge, $20 for the OEIS foundation in the honor of the solvers, I don't know the answer)
Discover an analog of Karatsuba's algorithm when you do
a*b=(a2*10^(2*d)+ a1*10^d+ a0)*(b2*10^(2*d)+ b1*10^d+ b0)=
(a2*b2)*10^(4*d)+ (a2*b1+a1*b2)*10^(3*d)+ ... + a0*b0
that naively needs nine multiplications at every recursive step, and see whether you can find
FIVE expressions of the form (Linear Combinations of a's)(Linear Combination of b's),
such that each of the 9 products a2*b2, ..., a0*b0,
can be written as linear combinations (with small integer coefficients) of these. In order to beat Karatsuba, you need
5 such "basis elements", since log[3](5)log[2](3) .
Added Jan. 27, 2014: See the students' nice
Solutions to
the 1st homework assignment
Programs done on Monday, Jan. 27, 2014
C2.txt ,
Contains procedures

B(f,x,a,b,eps): inputs an expression f in x, numbers a and b (a< b)
such that hopefully f(a) and f(b) have opposite signs,
and a small positive number eps, and uses the bisection method
to find an interval of length ≤ eps, where a root of f(x)=0 lies

NR(f,x,x0,N): performs N iterations in Newton's method to find an approximate
root of f(x)=0 starting with the initial guess x0.

TestNewton(f,x,x0,eps,N): inputs an expression
f in the variable x, an intial guess x0, a small eps
and a pos. integer N , and outputs the successive errors
abs(xnalpha) for n from 1 to N where alpha is the
appx. to the root of f(x)=0 found by the bisection method
Homework for Monday, Jan. 27, class (due Sunday Feb 2.,2014, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw2.txt. Please start the file with
#hw2.txt, Date, and YOUR NAME
If you do not want me to post your solutions, please write
#Do not post
Otherwise I will post them.

(For novices only)
Read and do all the examples, plus make up similar ones,
in pages 30 to 60 of Frank Garvan's
awesome Maple booklet
(
part 1,
part 2)
Only record a small sample in hw2.txt .
Note: a few commands are no longer valid in today's Maple. The
most important one is that " has been replaced by %.

(For everyone) Newton's method is the special case d=1 of
Householder's method, Adapt NR(f,x,x0,N) to write a procedure
Hd(f,d,x,x0,N)
that inputs a polynomial expression f in x, a positive integer d, a variable x,
an initial guess x0, and a positive integer N, and outputs the approximation
obtained by iterating N times procedure Hd1(f,d,x,x0); that would be the
analog of NR1
Hint: First evaluate the expression on the right side SYMBOLICALLY (as a rational function
of x_{n}), using Maple's bulitin differnetiation and algebraic manipulations
(namley the command "normal"). Also note that to differetiate an expression G of x,
d times, you do "diff(G,x$r);"

(For everyone) Adapt TestNewton to write TestHouseholder. Experiment with it
and see whether you get convergence rate of order d+1.

[Corrected Jan. 28, 2014, thanks to Andrew Lohr]
(for experts, optional for novices)
Write a program
Basins(c1,c2,c3,eps)
that inputs three integers c1, c2, c3 such that
c1 < c2 < c3
and a "resolution" eps, and outputs the list,
[S1,S2,S3]
such that if you apply Newton, say 5 times, (with floating point!)
to
f:=(xc1)*(xc2)*(xc3)
S1 is the set of numbers that are multiples of eps,
between c11 and c3+1 that converge to c1,
S2 is the set of numbers that are multiples of eps,
between c11 and c3+1 that converge to c2,
S3 is the set of numbers that are multiples of eps,
between c11 and c3+1 that converge to c3,
See if you can reproduce the data in
wikipedia article
about f(x)=(x+3)(x1)(x4).

(for everyone). Get ready for the next class by reading and understanding
Strassen's algorithm.
Added Feb. 3, 2014: See the students' nice
Solutions to
the 2nd homework assignment
Programs done on Thursday, Jan. 30, 2014
C3.txt ,
Contains procedures

Hil(n): the famous n by n Hilbert matrix

Per(M): inputs a square matrix (given as a list of lists)
and outputs the permanent of M (using the definition)

SD(M): inputs a square matrix (given as a list of lists)
and outputs the determinant using the definition
(VERY slow! do not use for matrices of dimension larger than 8)

sig1(p): the sign of a permutation p
((1) to the power the number of inversions of p)

AliceD(M): inputs a square matrix and outputs
the determinant using Dodgson (aka Lewis Carroll)
condensation method
Homework for Thursday, Jan. 30, class (due Sunday Feb 2.,2014, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw3.txt. Please start the file with
#hw3.txt, Date, and YOUR NAME
If you do not want me to post your solutions, please write
#Do not post
Otherwise I will post them.

(For novices only)
Read and do all the examples, plus make up similar ones,
in pages 60 to 90 of Frank Garvan's
awesome Maple booklet
(
part 1,
part 2)
Only record a small sample in hw3.txt .
Note: a few commands are no longer valid in today's Maple. The
most important one is that " has been replaced by %.

[For everyone]. Find out how dangerous floatingpoint calculations are by doing, say
for the Hilbert matrix Hil(50)
with(linalg): read `C3.txt`: Digits:=10: evalf(det(Hil(50))), det(evalf(Hil(50)));
Keep upping Digits:
with(linalg): read `C3.txt`: Digits:=20: evalf(det(Hil(50))), det(evalf(Hil(50)));
etc., until you get a good agreement. Record the smallest "Digits" that gives you good agreement.
Repeat it for Hil(n) for n=10,20,30,40, ..., 100, and record for each the smallest "good" Digits.
Can you predict a relation?

[For everyone] Modify procedure sig1(p) to write a program,inv1(p)
to count the number of inversions, and gf(n,q) to compute the
sum of q^{inv1(p)} over all the permutations of length n.
By factoring gf(n,q) for n=1..8, try to conjecture a "nice" formula for it.

[For experts, optional for novices] Write a procedure PerFast(M) that implements
Ryser's formula described in this
wikipedia article.
Compare the running time with Per(M) and PerFast(M) for random 8 by 8 and 9 by 9 matrices.

[Challenge, 2 dollars to the OEIS in honor of the people who would get it, no peeking!]
Let h(n)=det(Hil(n))
conjecture a formula for (h(n)/h(n1))/(h(n1)/h(n2)) for n ≥ 2, and deduce a closedform expression
(possibly involving factorials and producs thereof.)

[Challenge, 5 dollars to the OEIS in honor of the people who would get it, no peeking!]
Prove the formula rigorously, using Dodgson's rule.
Hint: first conjecture a more general frmula for the determinant of the generalized Hilbert matrix the n by nmatrix
given by
H[i,j]=1/(k+i+j1) \quad
and calling its determinant H(k,n), conjecture expressions for H(k,n)/H(k1,n) and H(k,n)/H(k,n1).
Then prove them by induction using Dogdson's condensation, and finally plug in k=0 to get h(n)
Added Feb. 3, 2014: See the students' nice
Solutions to
the 3rd homework assignment
Programs "done" (by Dr. Z.) on Monday, Feb. 3, 2014
Today was a "snow day", but in order that you would not get too sad, I wrote the code
that would have been written in class.
>C4.txt ,
Contains procedures

MatAdd(A,B): adds two matrices of the same size given as lists of lists

MatMult(A,B): multiplies two matrices A and B using the straightforward textbook definition

PadZeros(A): given a square matrix A pads it with zero to make its
dimension a power of 2

MatMultR(A,B): multiplies two square matrices A and B whose dimension is a power of 2
using (straightforward) recursion

Tra(A): the transpose of the matrix A

RandMat(m,n,R): a random m by n matrix with entries between R and R
given as a list of lists
Homework for Monday, Feb. 3, class (due Sunday Feb 9., 2014, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw4.txt. Please start the file with
#hw4.txt, Date, and YOUR NAME
If you do not want to have your homework posted please write
#Please do not post!

(For novices only)
Read and do all the examples, plus make up similar ones,
in pages 91 to the end of Frank Garvan's
awesome Maple booklet
(
part 1,
part 2)
Only record a small sample in hw4.txt .
Note: a few commands are no longer valid in today's Maple. The
most important one is that " has been replaced by %.

(For everyone)
Read and understand the maple code in
C4.txt .

(For everyone)
Extend file
C4.txt ,
and write a procedure MatMultRG(A,B), that uses MatMultR(A,B) in order to
multiply any two square matrices A and B of the same dimension, not
necessarily a power of 2.
[Hint, use PadZeros(A), and at the end of the day remove the spurious zeros]
Test it for random 5 by 5 matrices by comparing it with MatMul, and also with
the maple builtin matrix multiplication.

[Mandatory for experts, strongly recommended for novices]
Modify procedure MatMult(A,B) in
C4.txt
to write a procedure
Strassen(A,B) ,
that implements Strassen's famous (and beautiful!) algorithm
as it is described, e.g. in the
wikipedia article.
[obvious Hint: you may need to write a new procedure MatSubt(A,B) to
subtract matrices, or MinusMat(A) to multiply a matrix by 1].

Write a procedure RandMat(n,R) that inputs a positive integer n, and a large positive
integer R, and outputs a random n by n matrix with entries that are integers between 0 and R

Test your Srassen(A,B) vs. MatMultR(A,B) for random 64 by 64 matrices, 128 by 128, etc.
Do you notice any speedup?

[Optional]
So far no one did the challenge problem about conjecturing a closedform
for the determinant of the Hilbert matrix. Try again!
Added Feb. 10, 2014: See the students' nice
Solutions to
the 4th homework assignment
Programs done on Thursday, Feb. 6, 2014
C5.txt ,
Contains (in addition to those of C4.txt), procedures

RPow(A,k,x0): Raleigh quotient Power method for approximating
the dominant eigenvalue of a matrix A (given as a list of lists).
k is a positive integer indicating the number of guesses, and x0 is
an initial vector.

Ons1(z): the 2 by 2 transfer matrix for the 1 by infinity
Ising model columns are [1,1], rows are [1,1]
[Note this is correct, the largest eigenvalue is indeed z+1/z and this does not
have a 'singularily' (phase transition) unless z=0 (infinite temperature)
or z=inifnity (absolute zero), in other words water never boils in a onedimensional world]
Homework for Thurs., Feb. 6 class (due Sunday Feb 9., 2014, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw5.txt. Please start the file with
#hw5.txt, Date, and YOUR NAME
If you do not want to have your homework posted please write
#Please do not post!

By using RPow(A,k,x0) as a subroutine, write a procedure RPower(A,D1) that inputs a
square matrix, A, (given as a list of lists), and positive integer D1, and
finds an approximation for the dominant eignevalue with error that is believed to be less
than 10^{D1}, by doloping over k and quiting when the difference between
RPow(A,k,x0) and RPow(A,k+10,x0) is less than 10^{(D1+2)}, and
taking starting x0 either randomly, or "cleverly" (look it up).
Use floating point. So when you test it, use evalf(A) (it is much faster)

[For experts, optional to novices]
Write a procedure Onsm(m,z) that inputs a positive integer m and a number (or variable) z,
and outputs the famous 2^{m} by 2^{m} Onsager matrix defined as follows.
Both rows and columns are labeled by vectors of length m whose entries are {1,1}, where
the ith row (and column) corresponds to the {1,1} vector obtained by converting the integer i1
to binary (you can use base(i1,2);), and padding with 0's in front to make it of length m, and then replacing 0 by 1.
Then (if u and v are vectors as above)
M[u,v]:=z^(u[1]*v[1]+u[2]*v[2]+ ... + u[m]*v[m])*z^(u[1]*u[2]+u[2]*u[3]+ ... + u[m1]*u[m]+u[m]*u[1])

[For experts, optional to novices]
[Challenging, but you can do it!]
Let f_{m}(z) be the "logarithm(largest eigenvalue of Onsm(m,z))" divided by m.
write a program NumerOns(m,reso) that inputs a positive integer m, and outputs
the numerical values of f_{m}(z) for z between 0.2 and 3 in intervals of
reso (for exaple reso=0.1, or time permitting, reso=0.01). Confirm
Lars Onsager's seminal discovery
(see here)
of the existence of a "phase transition" by finding
(numerically!) the approximate locations of the maxima of f_{m}(z) as
m gets larger and larger. How big can you make m before the computer complains?
Added Feb. 10, 2014: See the students' nice
Solutions to
the 5th homework assignment
Special (paper!) Homework due Friday, February 14, 2014
Using xmaple, experiment, by trying out many numeric, a, c1, and c2, with the command
plots[implicitplot3d]( (x^2+y^2+z^21)^2a*(x^4+y^4)=0,x=c1..c1,y=c1..c1,z=c2..c2);
and find those parameters a,c1,c2 that makes it look at much as possible as a wrapped suckingcandy.
Then print it out, and write on top of it:
"Experimental Mathematics is as sweet as a candy"
and slide it under my door, Hill 704, on Friday, Feb. 14, 2014.
Programs done on Monday, Feb. 10, 2014
C6.txt ,
"The joy of guessing", contains procedures:

Hn(n): the determinant of the famous Hilbert n by n matrix (1/(i+j1))(1 <=i<=n)

HnS(N): the first N terms of Hn(n), starting at n=1

Hnr(n,r): the determinant of the generalized Hilbert n by n matrix (1/(i+j1+r))

HnrS(N): the first N terms of Hnr(n,r), starting at n=1

Ratio1(L): inputs a list L and outputs the list of L[i+1]/L[i]

GuessRF1(L,d,n): inputs a list of numbers (or expressions)
L, a nonneg. integer d, and a (discrete) variable name
n and tries to guess a rational function of n of degree
d, such that L[i]R(i)=0 for all i from 1 to nops(L)

GuessRF(L,n): inputs a list of numbers (or expressions)
L, and a (discrete) variable name
n and tries to guess a rational function of n of degree
≤ (nops(L)6)/2 such that L[i]R(i)=0 for all i from 1 to nops(L)
Homework for Monday, Feb. 10 class (due Sunday Feb 16., 2014, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw6.txt. Please start the file with
#hw6.txt, Date, and YOUR NAME
If you do not want to have your homework posted please write
#Please do not post!

let a(n) be a discrete function (e.g. a(n)=1/n in the case of the Hilbert matrix), define the
n by n matrix
H(n,r)[i,j]:=a(r+i+j1)) (1 ≤ i,j ≤ n)
and the double sequence
h(n,r)=det H(n,r)
Using Ratio1 (twice) and GuessRF, write a general procedure,
GuessRatioRatio(a,n,r)
that inputs the expression a, the discrete
variables n and r, and outputs a guessed expression, as a rational function of n, for
(h(n+2,r)/h(n+1,r))/(h(n+1,r)/h(n,r))
(in analogy with what we did in class).

Using the recurrence implied by Dodgson's rule
h(n,r)=(h(n1,r)*h(n1,r+2)h(n1,r+1)^2)/h(n2,r+2) ,
write a procedure ProveConj(Exp,n,r) that proves (rigorously!, by induction on n), that the output of
GuessRatioRatio(a,n,r) is indeed correct.

Try these out for a(n)=1/(n*(n+1)

The drawback of GuessRF(L,n) is that it can only handle rational functions where the degree of
the numerator is ≥ the degree of the denominator. By first writing
a procedure MyGuessRF1(L,d1,d2, n): that inputs, in addtion to the list L, pos. integers d1 and d2
and a variable n, and outputs a gnessed generating function with degree of top d1 and degree of
bottom d2, write a program MyGuessRG(L,n) that does not have the above drawback.

Write a program GuessDoubleRF(L,m,n) that inputs a list of lists L and variable names m,n, and
tries to find a rational function R in the variables m and n, such that for each
i,j, in the range L[i][j] equals R(i,j). Test it out for
L:=[seq([seq((i+j+1)/(i+j+4),j=1..10)],i=1..10)]
and make sure that you get
(m+n+1)/(m+n+4)
also try it out for other examples like that.
Added Feb. 17, 2014: See the students' nice
Solutions to
the 6th homework assignment
Programs "done" on Thurs., Feb. 13, 2014
Today was yetanothersnowday, so we must, once
again, resort to "distance learning". Only God can compute EXACTLY sin(x) for most x,
or exp(x), for most x, so we humans have to resort to
approximation. We already know from Calc2 that
the Taylor POLYNOMIALS at a certain place x=x0 do a pretty good
job for x NEAR x0, but what if you want to be fair, and approximate
the function, by a POLYNOMIAL (that we humans, and our computers can compute
exactly) UMIFORMLY on a specified interval? For this we have
POLYNOMIAL INTERPOLATION. The code that we would have
done today is in this file:
C7.txt ,
"Starting polynomial interpolation", that contains procedures:

InterPol(L,x):
inputs a list of pairs [[x1,y1],[x2,y2], ..., [xn,yn]] of numbers (or symbols)!
and a variable name x,and outputs the unique polynomial of degree ≤ n1,
let's call it P(x), such that
P(x1)=y1, P(x2)=y2, ..., P(xn)=yn,
For example, InterPol([1,3],x); should give 3.

InterPolF(f,x,L):
inputs an expression f, in the variable x, representing
a function, and a list L=[x1,...,xn], of numbers (or symbols)
outputs the unique polynomial of degree ≤ n1,
let's call it P(x), such that
P(x1)=f(x1), ..., P(xn)=f(xn). For example,
InterPolF(x^3,[1,2,3,4],x); should return x^3
Homework for Thurs., Feb. 13 class (due Sunday Feb 16., 2014, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw7.txt. Please start the file with
#hw7.txt, Date, and YOUR NAME

Carefully read, understand, and experiment with the two procedures,
InterPol(L,x), and InterPolF(f,x,L) of the file
C7.txt .

Scoop Joseph Lagrange, by running the purely routinelinearalgebra procedure
InterPol(L,x) on symbolic inputs, e.g.
f:=InterPol([[x1,y1],[x2,y2],[x3,y3]],x);
and then let Maple compute
factor(coeff(f,y1,1)); factor(coeff(f,y2,1)); factor(coeff(f,y3,1));
and analogously for
f:=InterPol([[x1,y1],[x2,y2],[x3,y3],[x4,y4]],x);
and CONJECTURE a GENERAL formula for the unique polynomial, P, of degree ≤ n1
such that P(xi)=yi

[Human] Rigorously prove this socalled Lagrange Interpolation Formula (no peeking!)
[Hint: By using the Euclidean algorithm for polynomials, and mathematical induction,
first prove a lemma: a polynomial of degree less than n that vanishes at n distinct places MUST be the zero polynomial]

Write a procedure
InterPolFast(L,x),
that inputs and outputs the SAME things as InterPol(L,x) but uses the readymade Lagrange Interpolation Formula
that you have just rediscovered. Compare the timing of

InterPolFast([seq([a[i],b[i]],i=1..10)],x); and InterPol([seq([a[i],b[i]],i=1..10)],x);

InterPolFast([seq([a[i],b[i]],i=1..20)],x); and InterPol([seq([a[i],b[i]],i=1..20)],x);

InterPolFast([seq([i,i^2],i=1..100)],x); and InterPol([seq([i,i^2],i=1..100)],x);

Realize the drawback of equidistance Lagrange interpolation by plotting (using Maple's plot):
n:=5:plot(eval(InterPolF(abs(x),x,[seq(i/n,i=n..n)]))eval(abs(x)),x=1..1);
and similarly for n=10, n=15, etc.

Read about the socalled Runge phenomenon in section 3.2.1. (pp. 5455) of
chapter 3 of the nice book written by Amparo Gil, Javier Segura, and Nico Temme,
and experiment with it, reproducing the graphs of figure 3.1 and do analogous things where the interval [5,5] is replaced by
[10,10] and [15,15].
Added Feb. 17, 2014: See the students' nice
Solutions to
the 7th homework assignment
Programs done on Monday, Feb. 17, 2014
C8.txt ,
The Lagrange Interpolation formula, NewtonCotes Quadrature, and starting Gaussian Quadrature.

LagP(L,x): inputs a list L=[x1,...,xn]
of numbers (or symbols)
and a variable name x, and outputs
the poly (xx1)*(xx2)* ...*x(xxn)

LIF(L,x): inputs a list of pairs:
L=[[x1,y1], ..., [xn,yn]]
and outputs
the unique polynomial in x of degree less than n
such that P(x1)=y1, ..., P(xn)=yn

FindNC(n): inputs a pos. integer n
and finds, using linear algebra,
the NewtonCotes
intermetiate points [x1,..., xn]
such that the integral of
f(x) from x=0 to x=1 equals
H[0]*f(0)+H[1]*f(1/n)+ ...+ H[n1]*f((n1)/n)+H[n]*f(1)
for all polynomials f(x) of degree ≤ n.

ApplyNC(f,x,n): inputs an expression f in the variable x,
and a positive integer n, applies the NewtonCotes formula
found by FindNC(n) to approximate the integral of f(x) from 0 to 1

FindGauss(n): inputs a pos. integer n
and finds, using nonlinear algebra
the Gauss quadrature formula
featuring a list of points
[x1,..., xn]
and a list of coefficients
[H[1], ..., H[n]]
such that
int(f(x),x=0..1)= H[1]*f(x1)+ ...+ H[n]*f(xn)
for ALL polynomials of degree less than 2*n.
It returns the pair of lists
[ [x1,..., xn], [H[1], ..., H[n]] ]
Homework for Monday Feb. 17 class (due Sunday Feb 23., 2014, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw8.txt. Please start the file with
#hw8.txt, Date, and YOUR NAME

Write a program
ApplyNCg(f,x,n,a,b) that uses NewtonCotes for approximating
int(f,x=a..b)
Hint: let Maple do a change of variable to the expression f to get it into
an integral from 0 to 1 and then use our ApplyNC(f,x,n)

Write a program
ApplyNCgg(f,x,n,a,b,K)
to approximate the integral of an expression f from x=a to x=b by chopping it into K equal intervals,
applying ApplyNCg to each subinterval, and adding them up.

Note that ApplyNCgg(f,x,n,a,b,K) uses roughly n*K function evaluations. Write a program
BestNC(f,x,M), that inputs an expression f of x, a positive integer M, and outputs
the pair of integers [n,K] such that M=n*K and
ApplyNCgg(f,x,n,0,1,K)
gets the closest to the true
int(f,x=0..1) ;
What are

BestNC(1/(1+x^2), x, 32);

BestNC(1/(1+x^4),x, 32);

BestNC(exp(x), x,32);

BestNC(log(1+x),x, 32);

Write a program
ApplyGaussian(f,x,n): inputs an expression f in the variable x, and a positive integer n,
and applies
FindGauss(n)
to approximate int(f,x=0..1);

Compare the errors of ApplyGaussian(f,x,n) and ApplyNC(f,x,n)
for
f=1/(1+x^2) f=1/(1+x^4) f=exp(x) f=log(1+x)
and n=1,2,3,4,5. Who is closer to the true value? By how much?
Added Feb. 24, 2014: See the students' nice
Solutions to
the 8th homework assignment
Programs done on Thursday, Feb. 20, 2014
C9.txt ,
Gaussian Quadrature the Gauss way (with a little help with Maple),
and naturally introducing orthogonal polynomials.
In addition to the relevant procedures from different classes, it
contains the following new procedures:

Pn(n,x): the unique momic polynomial of x of degree n
such that int(Pn(n,x)*x^i,x=0..1)=0 for i=0..n1

FindGaussSmart(n): Smart version of FindGauss(n)
inputs a pos. integer n
and finds, using nonlinear algebra
the Gaussian Quadrature numerical (appx.)
integration with clever
intermetiate points [x1,..., xn]
for int(f(x),x=0...1) that is EXACT
for polynomials of degree <=2*n1
int(f(x),x=0..1)= H[1]*f(x1)+ ...+ H[n]*f(xn). The output is the pair
[[x1,.., xn], [H[1], ..., H[n]]:

ApplyQ(f,x,Q)
inputs a function f of x, and a quadrature rule
Q=[X,H], computes the approximation to the integral that the quadrature rule is supposed to approximate.
Homework for Thurs. Feb. 20 class (due Sunday Feb 23., 2014, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw9.txt. Please start the file with
#hw9.txt, Date, and YOUR NAME
If you do NOT wish your homework to be posted please write
#Please do not post!
(or else I would post it).

[This problem is dedicated to my Patron Saint, GianCarlo Rota]
The mapping
P(x) > int(P(x),x=0..1) ;
is an example of LINEAR FUNCTIONAL (aka as Umbra)
on the (infinitedimensional) vector space of polynomials. By linearity, such
an umbra is completely specified by its action on monomials.
Write a program
ApplyUmbra(P,x,U,m)
that inputs a polynomial expression P in the variable x, a variable x, an expression U in m, and
a discrete variable m, where U describes the umbra
x^m > U(m)
and extended by linearity, and outputs the value when the umbra U is applied to P.
For example
ApplyUmbra(3+5*x+3*x^2,x,2^m,m)=3*2^0+5*2^1+3*2^2=25
ApplyUmbra(P,x,1/(m+1),m)= int(P,x=0..1)

Generalize procedure Pn(n,x) to write a procedure
PnG(n,x,U,m)
that inputs a nonneg. integer n, a variable name x, an expression U in the discrete variable m,
and the variable m, and outputs the unique monic polynomial,P, of degree n such that
when U is applied to the polynomials P*x^i (i=0..n1), you would get zero.
Check that PnG(n,x,1/(m+1),m); is the same as Pn(n,x) for n between 0 and 10.

By using GuessRF, write a procedure
GuessCoeff(n,i,U,m)
that inputs a discrete variable n, a nonnegative integer i, an expression U (describing an umbra)
in m, and a variable name m, and guesses a rational function in n for the
coefficient of x^(ni) in PnG(n,x,U,m).
What are

GuessCoeff(n,i,1/(m+1),m) for i=0..7

GuessCoeff(n,i,1/(m+3),m) for i=0..7

GuessCoeff(n,i,eval(m!),m) for i=0..7 ?

[Optional, a real challenge!] Write a procedure
GuessOP(n,r, U,m)
that inputs a SYMBOL (not a number) n, another discrete symbol r, an
umbra U (expression in m) and a variable name m, and outputs a closedform expression
in terms of the symbols n and r, for the coefficient of x^(nr) in PnG(n,x,U,m)
Added Feb. 24, 2014: See the students' nice
Solutions to
the 9th homework assignment
Programs done on Monday, Feb. 24, 2014
C10.txt ,
Orthogonal polynomials in general

AM(M,P,x): inputs a sequence of moments x^i>M[i+1]
outputs the "integral" of polynomial P of x

PnM(M,n,x)
inputs a sequence of numbers "the moments" such
that M[i+1] is what is mapped to x^{i}
outputs the nth degree MONIC polynomial in
the sequence of OPs such that
Umbra(P(n)*P(m))=0 if n &neq; m

PnMseq(M,n,x)

WtToM(w,n,x,a,b): the moments of f(x)>int(w(x)*f(x),x=a..b);

AnM(M,n): the first n terms of the A(n), featuring in the recurrence
x*P_{n}(x)= P_{n+1}(x)+A(n)*P_{n}(x)+B(n)*P_{n1}(x)
Reading Homework for Monday Feb. 24 class (due 11:50am, Thurs. Feb. 27, 2014)
Get ready for the next topic, numerical solutions of ordinary differential equations,
by reading sections 5.15.3 of Herb Wilf's
masterpiece,
Homework for Monday Feb. 24 class (due Sunday March. 2, 2014, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw10.txt. Please start the file with
#hw10.txt, Date, and YOUR NAME
If you do NOT wish your homework to be posted please write
#Please do not post!

[corrected and clarified, Feb. 25, 2014, thanks to Cole Franks]
Recall that procedure AnM(M,n) that we did in class, found the first n terms of the sequence
A(n)
featuring in the recurrence
x*P_{n}(x)= P_{n+1}(x)+A(n)*P_{n}(x)+B(n)*P_{n1}(x)
Wrire an analogous procedure, BnM(M,n) to compute the first n terms of the sequence B(n) above.

Using GuessRF, conjecture recurrences (i.e. EXPLICIT expressions for A(n), B(n) above) for

The Legendre polynomials, i.e.
M=[seq(1/(i+1),i=1..1000)]:

Laguerre polynomials:
M=[seq(i!,i=1..200)]:

M=[seq(1/(i+a),i=1..1000)]:
with (SYMBOLIC a)

Write a procdure (with option remember!)
PnFast(A,B,n,n1,x,a0)
That inputs expressions A, B in n, a positive integer n1, and a variable name x
and outputs the n1th term in the sequence of polynomials satisfying the recurrence
x*P_{n}(x)= P_{n+1}(x)+A(n)*P_{n}(x)+B(n)*P_{n1}(x)
under the initial conditions
P_{0}(x)=1 P_{1}(x)=x+a0
Hint: first rewrite the recurrence in a form expressing P_{n}(x) as a linear combination of
P_{n1}(x) and P_{n2}(x) .

Using PnFast(A,B,n,n1,x,a0) with the A and B that you found for the above moment sequences M (in problem number 2),
namely:

The Legendre polynomials, i.e.
M=[seq(1/(i+1),i=1..1000)]:

Laguerre polynomials:
M=[seq(i!,i=1..200)]:

M=[seq(1/(i+a),i=1..1000)]:
with (SYMBOLIC a)
find (time permitting!, the last one may take too long),
for each, the P_{1000}(x). Convince yourself that it would take for ever
to do it with PnM(M,1000,x) with the appropriate M.
Added March 3, 2014: See the students' nice
Solutions to
the 10th homework assignment
Programs done on Thurs., Feb. 27, 2014
C11.txt ,
Orthogonal polynomials in general

Orb(g,x0,K1,K2,d1): running the logistic map
x>g*x*(1x) K1+K2 times and displaying the last K2 terms
For example: Orb(2.5,.5,10000,1000);

OrbF(f,x,g,x0,K1,K2,d1):
running the general logistic map
x>g*f K1+K2 times and displaying
the last K2 terms. For example:
OrbF(2.5,.5,10000,1000);
masterpiece,
Homework for Thurs. Feb. 27 class (due Sunday March. 2, 2014, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw11.txt. Please start the file with
#hw11.txt, Date, and YOUR NAME
If you do NOT wish your homework to be posted please write
#Please do not post!

Do a second reading, and try to understand sections 5.15.3 of
Herb Wilf's
masterpiece.

By plain brute force, or by the "method of bisectionsearching", and using
nops({op(NumerBifOrbF(f,x,g,x0,K1,K2,d1))})
as an indicator of the probable orbitsize (it has to be a power of two), write a program
NumerBifOrbF(f,x,x0,K1,K2,d1,n,eps)
that inputs the same as OrbF(f,x,g,x0,K1,K2,d1),
(note, f(x) has to be a function that vanishes at 0 and 1)
as well as a positive integer n,
and a small positive number eps, and outputs a list, let's call it B, of length n,
such that B[i] is a small interval [ai,bi], with biai < eps, and such that
when g=ai the generalized logistic map has period 2^{i1} and
when g=bi it has period 2^{i} . For example,
NumerBifOrbF(x*(1x),x,.5,10000,1000,10,2,.02);
should yield something like : [[2.99,3.01], [3.44,3.46]] .
What are

NumerBifOrbF(x^2*(1x),x,.5,10000,1000,10,4,.02);

NumerBifOrbF(x*(1x)^2,x,.5,10000,1000,10,4,.02);

NumerBifOrbF(x^2*(1x)^2,x,.5,10000,1000,10,4,.02);

[Challenge, but please try!]
Write a program
FirstBif(f,x)
that inputs an expression f in x (that vanishes at x=0 and x=1), a varibale x, and a positive integer g
and finds the unique g1 such that
for g < g1 the map x>g*f(x) has (eventual) period one, and when g > g1 it has eventual period 2.
For example,
FirstBif(x*(1x),x);
should give 3.
What are

FirstBif(x^2*(1x),x);

FirstBif(x*(1x)^2,x);

FirstBif(x^2*(1x)^2,x);

FirstBif(x^3*(1x)^3,x);

Look up the definition of the Mandelbrot set, and write a program
PlotMS(reso,K);
that inputs a small resolution reso, and a positive integer K, and includes those c
in the complex plane (approximated with a grid of resolution reso)
such that iterating
z>z^2+c
(starting at z=0), K times does NOT go far away.
Save the output of
PlotMS(.05,10);
as .jpg file and put it in your directory, as MS.jpg
[Hint: look up the plot commands for plotting a set of points, and for saving as a .jpg file]
Added March 3, 2014: See the students' nice
Solutions to
the 11th homework assignment
Programs done on Monday, March 3, 2014
C12.txt ,
Elucidating Picard's proof of the existence and uniquness of an initial value problem
for an ordinary differential equation, and trying out a powerseries method for
approximating solutions, in view for their numerical solutions.
Added March 7, 2014: Andrew Lohr pointed out two shorcomings of procedure
Z1 (and hence Z). These have been corrected. Please discard old version.
Contains the following procedures:

DS(f,x,y,x0,y0): inputs an expression, f in the variables x and y, uses Maple's dsolve
to solve the initial value problem
y'(x)=f(x,y(x)), y(x0)=y0
If Maple can't do (usually the case), it returns FAIL.

Picard(f,x,y,x0,y0,N): the first N+1 terms
in the Picard sequence that is supposed to
converge to the solution of y'(x)=f(x,y(x)), y(x0)=y0
Warning: takes a long time for large, N, seems to be only of theoretical value, in general.

One step in the "Zeilberger" naive approach, using power series to find
successive MacLaurin polynomials of solutions of initial value problem diff.eqs. of the form
y'(x)=f(x,y(x)), y(0)=y0 .

Z(f,x,y,y0,N): the list of the first N+1 MacLaurin polynomials x of sol. of
y'(x)=f(x,y(x)), y(0)=y0

ERC(f,x): The empirical radius of convergence, using the ratio test, of the Maclaurin polynomial f in x
(f has to be of at least degree 30 for this to be useful. It also returns the
approximation due to the MacLaurin polynomial of degree one less.
If the two values are close, it is probably a good estimate.
Reading Homework for Monday March 3, 2014, class (due 11:50am, Thurs. March 6, 2014)
Read section 5.5 of Herb Wilf's
masterpiece,
Homework for Monday March 3, class (due Sunday March. 9, 2014, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw12.txt. Please start the file with
#hw12.txt, Date, and YOUR NAME
If you do NOT wish your homework to be posted please write
#Please do not post!

Modify procedure
Z(f,x,y,y0,N)
[Note new version, thanks to Andrew Lohr!]
calling it
Zg(f,x,y,y0,N)
so that f can be ANY analytic function of x and y (i.e. one that posseses a Taylor series
around the origin (0,0). (Hint, use the "taylor" command (in one variable , x))
Note: if f is a polynomial of x and y, then it should give the SAME output (test it first to make sure that you get the same output)
What are:

Zg(cos(x*y)+exp(x),x,y,1,10)[11];

Zg(exp(sin(x)+y)+1+x+y,x,y,1,10)[11];

Zg(cos(x*y)+exp(x)+1,x,y,1,10)[11];

Modify Zg(f,x,y,y0,N) to write a procedure
PolAppx(f,x,y,y0,x1,N,eps):
that inputs

(i)an expression f in x and y that is analytic in x and y near the origin (0,0)

(ii) variable names x,y

(iii) a number y0

(iv) a number x1

(v) an integer N

a small positive number eps
and first checks that abs(x1) is less than the ERC(Zg(f,x,y,y0,N)/2 (to be safe!),
and if does not, returns FAIL. Then finds the MacLaurin polynomial of least degree
such that
subs(x=x1,Zg(f,x,y,x0,N)[N+1]  Zg(f,x,y,x0,N)[N]) < eps/100
followed by its value at x=x1.
What are:

PolAppx(1+x^2+y^3,x,y,1,0.3,50,10^(10));

PolAppx(cos(x)*exp(y)+1+x+y,x,y,1,0.3,50,10^(10));

[Human] Prove that if f is a polynomial in x and y with INTEGER coefficients, then
the the Maclaurin series of the solution of
y'(x)=f(x,y(x)) , y(0)=1
can be written as
Sum(c(n)*x^n/n!,n=0..infinity)
where c(n) are integers!

[Challenge!] Modify procedure Z(f,x,y,y0,N) (and hence Z1) to write a procdure
Zm(f,x,y,y0,N)
where now the input is

(i) a LIST f of polynomials in x and the variables y[1],y[2], ..., y[nops(y)] (the dependent variable), where y is a list of variables (see below)

(ii) x a singlevariable (the independent variable)

(iii) y a LIST of depenendent variables (of the same length as the list f)

(iv) y0 is a LIST of numbers

(v) N is a positive integer
and outputs the LIST that is an approximate solution (up to degree N) of the system of diff.eqs.
y[1]'(x)=f[1](x,y[1], y[2], ..., y[nops(y)])
....
y[nops(y)]'(x)=f[nops(f)](x,y[1], y[2], ..., y[nops(y)])
Subject to the initial conditions
y[1](0)=y0[1] , ...., y[nops(y)](0)=y0[nops(y)] , ....,
In particular, if f is a polynomial in x,y, and y is a variable, and y0 is a number
Zm([f],x,[y],[y0],N);
should be the list of length n, [POL], where POL is the output of Z(f.x.y,y0,N)[N+1];
What are:

Zm([1+x*y1+y2, 3+x*y1^2+y2],x,[y1,y2],[1,2],10);

Zm([x+y1,y1*y2+x,x+y1*y2*y3],x,[y1,y2,y3],[1,2,1],10);
Added March 10, 2014: See the students' nice
Solutions to
the 12th homework assignment
Programs done on Thurs., March 6, 2014
In preparation for Pi Day we did
C13.txt ,
that contains the following procedures:

JesusPi1(n): computes the nth partial sum of the first formula (for 128/Pi^{2})
in Jesus Guillera's beautiful homepage.

JesusPi(n): computes an approx. to Pi from the above appx. to 128/Pi^2
We also started numerical solutions of odes:
C13a.txt ,
that contains procedures:

Euler1(f,x,y,x0,y0,h): One step in Euler's method to
approximate solutions of y'(x)=f(x,y(x)), y(x0)=y0 with
mesh size h

Euler(f,x,y,x0,y0,x1,h): Euler's method approximating
y'(x)=f(x,y), y(x0)=y0, all the way to x1
with increments (mesh size) h

impEuler1(f,x,y,x0,y0,h): One step in the improved Euler's method to
approximate solutions of y'(x)=f(x,y(x)), y(x0)=y0 with
mesh size h

impEuler(f,x,y,x0,y0,x1,h): The Improved Euler's method approximating
y'(x)=f(x,y), y(x0)=y0, all the way to x1
with increments of h
Homework for Thurs. March 6, class (due Sunday March. 9, 2014, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw13.txt. Please start the file with
#hw13.txt, Date, and YOUR NAME
If you do NOT wish your homework to be posted please write
#Please do not post!

Write procedure OtherJesusPi(n) implementing the second (still conjectured)
infinite series for Pi in
Jesus Guillera's beautiful homepage
for 128*sqrt(5)/Pi^2 . Check that indeed every new terms gives 5 more digit
Hint: do not compute (6*n)!/n!^6 directly, but rather use an analogous approach to the one we did in class,
that updates the current summand in terms of the previous one.

Write a procedure
RK4(f,x,y,x0,y0,x1,h)
Implementing the RungeKunge Method (analogously to Euler and Improved Euler) as it
is described, for example in
this handout

Write three procedures,
EulerH(f,x,y,x0,y0,x1,Di); impEulerH(f,x,y,x0,y0,x1,Di); and RK4H(f,x,y,x0,y0,x1,Di);
that inputs an expression f in x, y such that Maple's dsolve knows to find the
exact solution to theinitial value problem
y'(x)=f(x,y(x)) , y(x0)=y0 ,
(if it can't the procedure should return FAIL), and then
finds out the largest mesh size h for which the Euler, Improved Euler, and RungeKutta (RK4)
gives you an approxination correct up to Di digits.
What are

EulerH(y,x,y,0,1,1,10); , impEulerH(y,x,y,0,1,1,10); , RK4H(y,x,y,0,1,x1,10);

EulerH(1+x+y,x,y,0,1,1,10); , impEulerH(1+x+y,x,y,0,1,1,10); , RK4H(1+x+y,x,y,0,1,x1,10);

EulerH(1+x*y,x,y,0,1,1,10); , impEulerH(1+x*y,x,y,0,1,1,10); , RK4H(1+x*y,x,y,0,1,x1,10);

Do a first reading of
the Wikipedia entry on RungeKutta methods.
Added March 10, 2014: See the students' nice
Solutions to
the 13th homework assignment
Programs done on Monday, March 10, 2014
C14.txt :
General RungeKutta methods, following
wikipedia.
Note: two new procedures, GRK1s and ProveB were added after class
was over, please read and understand them. GRK1s is a symbolic
adaption of GRK1, where the input h is symbolic, and in order
to avoid Maple crashing (or taking for ever), it truncates
the intermediate polynomials in h, up to the specified order,
and that's all you need for the required proof, accomplished
in the new procedure ProveB.

RK4: NOT a procedure, but just an object,
the Butcher representation of the classical RungeKutta
[[],[1/2],[0,1/2],[0,0,1]],[1/6,1/3,1/3,1/6],[0,1/2,1/2,1]:

GRK1(a,b,c,f,x,y,x0,y0,h): One step in the
General RungeKutta with Butcher matrix a,b,c
(see wiki)

GRK(a,b,c,f,x,y,x0,y0,h,x1): The general RungeKutta
with Butcher matrix a,b,c appx. to
the ode
y'(x)=f(x,y(x)), y(x0)=y0
all the way to x1 with stepsize h
the output is pair of lists
Xvalues, Yavlues

[From C12.txt]:
Z(f,x,y,y0,N): the Nth MacLaurin polynomial
in x of the sol. of y'(x)=f(x,y(x)), y(0)=y0

P(x,y,a,d): a generic polynomial in x,y of degree d
in variables x,y and degree d, with symbolic coefficients of the form a[i,j],
where a is a letter.

[Added after class was over]
GRK1s(a,b,c,f,x,y,x0,y0,h,R): Symbolic rendition (with symbolic h),
of one step in the
General RungeKutta with Butcher matrix a,b,c,
BUT, only retaining powers of h up to R. Used for
proving RIGOROUSLY, that the given RungeKutta method a,b,c
is of order R. (see below)

ProveB(B,d): inputs a Butcher triplet [a,b,c]
and a positive integer d, rigorously decides whether the
explicit RungeKutta method, based on B is of order d. Try:
ProveB([ [[],[1/2],[0,1/2],[0,0,1]],[1/6,1/3,1/3,1/6],[0,1/2,1/2,1] ],4);
Note:Now the Butcher matrix, B, is given as a list [a,b,c], rather then,
as simply a,b,c like in GRK1 and GRK1s.
Homework for Monday March 10, class (due Wed. March. 12, 2014, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw14.txt. Please start the file with
#hw14.txt, Date, and YOUR NAME
If you do NOT wish your homework to be posted please write
#Please do not post!
NOTE DUE DATE, It is March 12, 2014, SO THAT YOU
WOULD HAVE A RELAXING SPRING BREAK! Of course, you are welcome to finish it
earlier!

Carefully study the two procedures,
GRK1s(a,b,c,f,x,y,x0,y0,h,R)
and
ProveB(B,d)
and understand them. Convince yourself, that if ProveB(B,d) returns true,
then you have a RIGOROUS proof that the method given by the Butcher triplet,
is of order d.

Check that indeed
ProveB([ [[],[1/2],[0,1/2],[0,0,1]],[1/6,1/3,1/3,1/6],[0,1/2,1/2,1] ],4);
gives you true. Then modify the weights
[1/6,1/3,1/3,1/6]
of the weightedaverage,e.g.
[1/4,1/4,1/4,1/4] , [1/10,2/10,3/10,4/10] , [1/10,4/10,4/10,1/10]
,
(keeping everything else the same), etc., and see if any of them has order 4,
and if not, what is the correct order (0,1,2, or 3).

Write a procedure, DelinZheng(lambda), that inputs a number lambda,
and outputs the RungeKutta method described at the bottom of p.1 and the
top of p.2 of the
wikipedia article,
as a triplet [a,b,c], the socalled Butcher matrix. Make sure
that DelinZheng(2)=[RK4].

Using procedure ProveB, prove rigorously that with lambda=1,3,4,5 you also
get fourthorder methods. What are the orders for other lambda?
For example, lambda=2.5, 3.5?

To prove that a given RungeKutta method is definitely of order d, you need
to input a generic polynomial f of degree d in x,y,
with symbolic coefficients,
but to disprove it, any numeric such polymial would do.
Also to semirigorously prove it (much faster, for larger d), you can
test it for numerous random such polynomials. In view of that,
modify procedure
P(x,y,a,d)
and write a procedure
RP(x,y,M,d)
that inputs variables x,y, a positive integer M, and a positive integer d
and outputs a polynomial in x,y, of (total) degree with numeric, integer coefficients
between 1 and M
Hint: use rand(1..M)()

Write a procedure
srProveB(B,d,M,K)
that rigorously proves that B is NOT of order d if it returns false,
and semirigorously proves that B is of order d if it returns true,
by checking K different randomly chosen polynomials
generated by RP(x,y,M,d).
Added March 13, 2014: See the students' nice
Solutions to
the 14th homework assignment
Programs done on Thurs., March 13, 2014
C15.txt:
Starting to discover RungeKutta methods ab initio.
[The rest would be done by you in the homework].
In addition the old procedures from C14.txt, it contains the following procedures.
Homework for Thurs. March 13, class (due Sun. March. 16, 2014, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw15.txt. Please start the file with
#hw15.txt, Date, and YOUR NAME
If you do NOT wish your homework to be posted please write
#Please do not post!

Write a program FindRK(s,d) that inputs positive integers s and d and finds all
RungeKutta methods with s steps and order d. For example,
RK(4,4)
should either give DelinZheng, or something even more general (i.e. with more than one parameter)
for which DelinZheng would be a special case.
Warning: I have not tried it, a priori you may get a system of algebaric equations that may be
too hard for Maple, but you should try it, and if it crashes, or takes too long, try to do things more
cleverly/modestly. For example, you may want to have the clist part of the input, and
try out may c's, for example [0,1/2,1/2,1].
Also taking random polynomials for f may be too cumbersome. Why not take as simple f's as possible
for getting a system of algebraic equations for the unknown entires of the Butcher triplet [a,b,c].
Also, start pennypinching, and use the fact that the sum of the b's should be 1.

[Challenge, I don't know the answer] Without peeking in the internet, find an order5 RungeKutta method,
using the present approach, possibly enhanced.
[Added March 15, 2014, 4:38pm. According to L. Nick Trefethen's essay in
p.604615 of the "Princeton Companion of Mathematics" (W.T. Gowers, ed.), available from this
secret url,
there is no order5 RungeKutta with 5 stages (i.e. s=5). Can you have Maple prove it rigorously?
Can you find (w/o peeking) an order5 RK with s=6?
[It may be a heavy computation, but maybe with cleverness you can do it!]
Added March 17, 2014: See the students' nice
Solutions to
the 15th homework assignment
Programs done on Monday, March 24, 2014
C16.txt:
Running and discovering Multistep Adamsbashford methods; Implicit RungeKutte methods.
Contains the following procedures (and Butcher matrices RK4, GL)

The Butcher matrix of (the explicit) RK4:
RK4:= [[0,0,0,0],[0,0,0,1/2],[0,0,0,1/2],[0,0,0,1]],[1/6,1/3,1/3,1/6],[0,1/2,1/2,1]:

The Butcher matrix of (the implicit) GaussLegendre RungeKutta method:
GL:=[[1/4,1/41/6*sqrt(3)],[1/4+1/6*sqrt(3),1/4]],[1/2,1/2],[1/2sqrt(3)/6,1/2+sqrt(3)/6]:

AAB(f,x,y,x0,y0,h,x1,L), inputs

An expression f in variables x, and y
(describing the function (x,y)>f(x,y)

initial conditions y(x0)=y0

h: mesh size

x1 last point needed

L a list of numbers describing the AdamsBashford method
y(x+h) ~ y(x)+ h*add(L[i+1]*f(xh*i,y[xh*i]),i=0..nops(L)1)
(e.g. Euler is L=[1], 2step AdamsB: L:=[3/2,1/2])
Outputs approximations to the solutions of the initial value ODE
y'(x)=f(x,y(x)), y(x0)=y0
starting at x0, and ending at x1, with increments of h. It returns the list of the x's, followed by the list of the approximate y's.

FAB(d): discovers the AdamsBashdford dstep method (of order d)

GIRK1(a,b,c,f,x,y,x0,y0,h): One step at the Implicit
RungeKutte method with Butcher matrix a,b,c

GIRK(a,b,c,f,x,y,x0,y0,h,x1): The general RungeKutta
with Butcher matrix a,b,c appx. to
the ode
y'(x)=f(x,y(x)), y(x0)=y0 all the way
to x1 with stepsize h
the output is pair of lists
Xvalues, Yvalues
Note that we now need a different convention for the matrix a, it is a list of length s whose entries are s by s matrix.
So instead of the RK4 defined in C14.txt, we have this one:
RK4:= [[0,0,0,0],[0,0,0,1/2],[0,0,0,1/2],[0,0,0,1]],[1/6,1/3,1/3,1/6],[0,1/2,1/2,1]:
given in the current file C16.txt
Reading Homework for Monday, March 24, class (due Thurs. March. 27, 2014, 11:59am)
Read and understand
Programming Homework for Monday, March 24, class (due Sun. March. 30, 2014, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw16.txt. Please start the file with
#hw16.txt, Date, and YOUR NAME
If you do NOT wish your homework to be posted please write
#Please do not post!

Write a procedure
BestAdams(f,x,y,x1)
that inputs an expression f in x, and y such that
dsolve({diff(y(x),x)=f,y(0)=1},y(x));
gives you a closedform answer, and outputs the d in {1,2, ..., 10} and h in {1/10,1/100,1/1000}
such that
AAB(f,x,y,0,1,h,x1,FAB(d))
is as close as possible to the true y(x1), and also returns the error.
What are

BestAdams(y,x,y,1)

BestAdams(x+y,x,y,5)

BestAdams(x+3*y,x,y,10)

Even though
RK4:= [[0,0,0,0],[0,0,0,1/2],[0,0,0,1/2],[0,0,0,1]],[1/6,1/3,1/3,1/6],[0,1/2,1/2,1]:
is order 4 while
GL:=[[1/4,1/41/6*sqrt(3)],[1/4+1/6*sqrt(3),1/4]],[1/2,1/2],[1/2sqrt(3)/6,1/2+sqrt(3)/6]:
is of order 2, the latter may give better results. Find examples where
GL gives you much better approximations than RK4, due to its superior stability.
Added March 31, 2014: See the students' nice
Solutions to
the 16th homework assignment
Programs done on Thurs., March 27, 2014
C17.txt:
getting ready for numerical solutions of PDEs, but for now
staying in one dimension, as a warmup.

GPi(N): the geometrical discrete Pi with meshsize N,
i.e. the number of lattice points (x,y) such
that x^2+y^2
Dsin(x,h): the discrete analog of the sine function with meshsize h

Srec(L,Ini,n): The nth term of the
solution of the linear recurrence equation
with constant coefficients
a(n)=L[1]*a(n1)+L[2]*a(n2)+ ... + L[1]*a(nnops(L))
and
a(0)=Ini[1], ..., a(nops(Ini)1)=Ini[1]

dBVP(L,Ini,Fini,N): dBVP: Discrete boundary value problem
The list of length N+1
giving the values of a(n) for n=0 to n=N
to the solution of the linear recurrence equation
with constant coefficients
a(n)=L[1]*a(n1)+L[2]*a(n2)+ ... + L[1]*a(nnops(L))
a(0)=Ini[1], ..., a(nops(Ini)1)=Ini[1].
a(N)=Fini[1], a(N1)=Fini[2], ..
Homework for Thurs., March 27, class (due Sun. March. 30, 2014, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw17.txt. Please start the file with
#hw17.txt, Date, and YOUR NAME
If you do NOT wish your homework to be posted please write
#Please do not post!

[Corrected March 29, 2014, 6:40pm, thanks to Anthony Zaleski]
Extend procedure
Srec(L,Ini,n):
and write procdure
SIrec(L,Ini,f,m,n)
with the same input as Srec but in addition an expression, f, in the
discerte variable m,
that solves Inhomogeneous linear recurrences
a(n)=L[1]*a(n1)+L[2]*a(n2)+ ... + L[1]*a(nnops(L))+f(n)

Do the analogous thing to
dBVP(L,Ini,Fini,N):
writing a procedure
dIBVP(L,Ini,Fini,f,n,N):

Use procedure dBVP(L,Ini,Fini,N)
to write a procedure
GR(p,N)
that inputs a real number p, between 0 and 1 (or symbolic), and
a positive integer N, and outputting a list of length N+1, let's call
it G, such G[i+1] is your chance of exiting a winner in
a casino with the following rule:
if you enter a casino
[Added March 30, 2014, thanks to Anthony Zaleski], With i dollars
where at each step you win
a dollar with probability p and lose a dollar with probability 1p, and
you stop playing when you are either broke (0 dollars), or got N dollars
[Hint: first do the human part, find a recurrence for that probability,
taking into account that if currently you have i dollars, then next time
you have i+1 dollars with probability p and i1 dollars with probability 1p.

Use the procedure that you wrote,
dIBVP(L,Ini,Fini,f,n,N)
to write a procedure
LGR(p,N)
that inputs the same thing as above , and outputs the list whose ith entry is the EXPECTED duration of
the game, if you currently have i dollars, until you end it (either winner (with N dollars) or loser with 0 dollars).
[Hint: first do the human part, find an inhomogeneous recurrence for that expected
duration, using conditional probability,
taking into account that if currently you have i dollars, then next time
you have i+1 dollars with probability p and i1 dollars with probability 1p,
but you have wasted one round.

Using the outputs of GR(1/2,N), and LGR(1/2,n), conjecture closedform
expressions for the probability of winning, and expected duration, respectively,
if you currently have i dollars, and you are in a fair casino (of course, this
is fictional! a fair casino is an oxymoron), hence p=1/2, and you have to exist
as soon as you are broke, or have N dollars.

Do first reading of chapter 3 of
chapters 35 of Garrett Birkhoff's classic.
Added March 31, 2014: See the students' nice
Solutions to
the 17th homework assignment
Programs done on Monday, March 31, 2014
C18.txt:
getting ready for numerical solutions of PDEs, but
still staying in one dimension, as a warmup.

SrecS(L,Ini,N): inputs a linear recurrence equation
with constant coefficients, coded by L
and a list of initial values, Ini
outputs the LIST of length N+1 such that
L[i+1]=Srec(L,Ini,n)

dBVPc(L,Ini,Fini,N): a (hopefully) clever way
to dBVP(L,Ini,Fini,N). It outputs
he list of length N+1
giving the values of a(n) for n=0 to n=N
to the solution of the linear recurrence equation
with constant coefficients
a(n)=L[1]*a(n1)+L[2]*a(n2)+ ... + L[1]*a(nnops(L))
a(0)=Ini[1], ..., a(nops(Ini)1)=Ini[1]
a(N)=Fini[1], a(N1)=Fini[2], ..

P1ch(f,x,h,mu0,mu1): the exact solution of the 1D Poisson eq.
u''(x)=f(x), u(0)=mu0, u(1)=mu1
where f is an expression in the variable x,
at 0,h,2*h, ..., 1

P1d(f,x,h,mu0,mu1): approximate solution of the 1D Poisson eq.
u''(x)=f(x), u(0)=mu0, u(1)=mu1
where f is an expression in the variable x,
h the mesh size, and using
u''(x) ~ u(x+h)2*u(x)+u(xh)

GlE(f,x,h,mu0,mu1):
Global error in using a mesh size h
to solve the 1D Poisson eq.
u''(x)=f(x), u(0)=mu0, u(1)=mu1
Homework for Monday, March 31, class (due Sun. April 6, 2014, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw18.txt. Please start the file with
#hw18.txt, Date, and YOUR NAME
If you do NOT wish your homework to be posted please write
#Please do not post!

Using SIrec(L,Ini,f,m,n) that you wrote for hw17.txt,
(see e.g., Nathan Fox's homework 17),
write a oneline procedure
SIrecS(L,Ini,f,m,N)
analogous to SrecS(L,Ini,N).

Using the abovementioned SIrecS(L,Ini,f,m,N), that you just wrote, adapt
dBVPc(L,Ini,Fini,N)
that we did in class (see C18.txt)
to write a procedure
dBVPIc(L,Ini,Fini,f,m,N)
for cleverly solving boundaryvalue recurrence equations.

Using the above procedure, write a procedure
P1dClever(f,x,h,mu0,mu1)
that does the same thing as procedure P1d(f,x,h,mu0,mu1) written in class, but hopefully much faster.

Adapt GlE(f,x,h,mu0,mu1) to write
GlEclever(f,x,h,mu0,mu1)
that takes advantage of P1dClever(f,x,h,mu0,mu1)

Write a procedure
EmpiricalStability(f,x,mu0,mu1, ListH)
that finds the equence of ratios
GlEclever(f,x,h,mu0,mu1)/h^2
For all h in ListH.
What are

EmpiricalStability(x^4,x,0,1,[1/10,1/20,1/30,1/40]);

EmpiricalStability(exp(x),x,0,1,[1/10,1/20,1/30,1/40]);

EmpiricalStability(1/(3x),x,0,1,[1/10,1/20,1/30,1/40]);
Added April 7, 2014: See the students' nice
Solutions to
the 18th homework assignment
Programs done on Thurs., April 3, 2014
C19.txt:
Numerical solutions of the 1Dimensional Heat Equation, following the
wikipedia article of Finite Difference Methods
Contains procedures:

Heq1(U,x,h,k,x1,t1): The Explicit method
inputs an expression U of the variable x
two small numbers h and k (the mesh sizes in x and t resp.)
and numbers x1, t1, where x1 is between 0 and 1
and a final time t1, an
outputs an approximation to the value of u(x1,t1),
where
u(x,t) is the solution of the 1DHeat equation
u_{xx}=u_t, u(x,0)=U(x), u(0,t)=0, u(1,t)=0

Heq1f(U,x,h,k,x1,t1): Floatingpoint version of the above

Heq2(U,x,h,k,x1,t1): same as Heq1(U,x,h,k,x1,t1)
but using the implicit method.
Homework for Thurs., April 3, class (due Sun. April 6, 2014, 10:00pm)
All homework should be uploaded to the secret directory that you gave me,
as file hw19.txt. Please start the file with
#hw19.txt, Date, and YOUR NAME
If you do NOT wish your homework to be posted please write
#Please do not post!

Adapt procedure
Heq2(U,x,h,k,x1,t1)
to write procedure
Heq3(U,x,h,k,x1,t1)
that implements the CrankNicolson method descibed in the abovementioned
wikipedia article

Write a procedure
Compare(U,t,x,h,k,x1,t1)
that inputs an expression U in the variable t and x, that you know beforehand
satisfies the Heat Equation and the above boundary and initial conditions
(it should return FAIL if it does not!). For example any
finite linear combination of
sin(n*Pi*x) * exp(n^2*Pi^2*t)
[typo corrected, thanks to Anthony Zaleski]
for integer n, and ouputs, a list with FOUR values
[The Exact Value, OutputOfHeEq1, OutputOfHeEq2,OutputOfHeEq3]
and an integer in the set {1,2,3} indicating who got closest to the exact value.

[Challenge] By approximating the Laplace operator
diff(u,x,x)+ diff(u,y,y)
by the socalled fivepoint stencil
(u(x+h,y)2*u(x,y)+u(xh,y))/h^2 + (u(x,y+k)2*u(x,y)+u(x,yk))/k^2
write a procedure
Dirichlet(f,x,y,h,k,x1,y1)
that inputs an expression f in the variables x and y, small positive numbers h and k, and number x1, y1 between 0 and 1
(that must be multiples of h and k respectively), and outputs a finitedifference approximation to
u(x1,y1), using mesh size h in the x direction, and meshsize k in the ydirection, and u(x,y) is
the solution of the Poisson equation
diff(u,x,x)+ diff(u,y,y)= f
in the unit square [0,1]x[0,1], that vanishes on the four sides of the unit square.
What are

Dirichlet(x^2+y^2,x,y,1/10,1/10,1/2,1/2)

Dirichlet(x^2+y^2,x,y,1/20,1/20,1/2,1/2)

Dirichlet(x^2+y^2,x,y,1/30,1/30,1/2,1/2)
Added April 7, 2014: See the students' nice
Solutions to
the 19th homework assignment
Programs done on Mon., April 7, 2014