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


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 computer-algebra 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 (Runge-Kutta) 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 computer-generated and computer-assisted 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. 604-615 from this url.

How to submit your homework

If you don't already have a public_html directory, go to command-line Unix (or Linux), and type

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

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
  1. (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 %.

  2. (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.

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

  4. (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^2-2,x,1,10^(-10));
    should be the same as evalf(sqrt(2));

  5. (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 i-th entry is the average running time of M(a,b) (and K(a,b)) respectively of ten randomly-drawn pairs of integers with d digits. See whether they are asymptotically a constant times d2 and d(log[2](3))= d1.584962501

  6. (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

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.
  1. (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 %.

  2. (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 xn), using Maple's bulit-in 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);"

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


  4. [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:=(x-c1)*(x-c2)*(x-c3)
    S1 is the set of numbers that are multiples of eps, between c1-1 and c3+1 that converge to c1,
    S2 is the set of numbers that are multiples of eps, between c1-1 and c3+1 that converge to c2,
    S3 is the set of numbers that are multiples of eps, between c1-1 and c3+1 that converge to c3,

    See if you can reproduce the data in wikipedia article about f(x)=(x+3)(x-1)(x-4).

  5. (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

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.
  1. (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 %.

  2. [For everyone]. Find out how dangerous floating-point 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?

  3. [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 qinv1(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.

  4. [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.

  5. [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(n-1))/(h(n-1)/h(n-2)) for n ≥ 2, and deduce a closed-form expression (possibly involving factorials and producs thereof.)

  6. [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 n-matrix given by
    H[i,j]=1/(k+i+j-1) \quad
    and calling its determinant H(k,n), conjecture expressions for H(k,n)/H(k-1,n) and H(k,n)/H(k,n-1). 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

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!

  1. (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 %.

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

  3. (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 built-in matrix multiplication.

  4. [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].

  5. 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

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

  7. [Optional] So far no one did the challenge problem about conjecturing a closed-form 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

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!

  1. 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 do-loping 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)

  2. [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 2m by 2m Onsager matrix defined as follows.

    Both rows and columns are labeled by vectors of length m whose entries are {-1,1}, where the i-th row (and column) corresponds to the {-1,1} vector obtained by converting the integer i-1 to binary (you can use base(i-1,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[m-1]*u[m]+u[m]*u[1])

  3. [For experts, optional to novices]
    [Challenging, but you can do it!]
    Let fm(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 fm(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 fm(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^2-1)^2-a*(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 sucking-candy. 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:

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!

  1. 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+j-1)) (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).

  2. Using the recurrence implied by Dodgson's rule
    h(n,r)=(h(n-1,r)*h(n-1,r+2)-h(n-1,r+1)^2)/h(n-2,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.

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

  4. 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.

  5. 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 yet-another-snow-day, 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:

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

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

  2. Scoop Joseph Lagrange, by running the purely routine-linear-algebra 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 ≤ n-1 such that P(xi)=yi

  3. [Human] Rigorously prove this so-called 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]

  4. Write a procedure

    InterPolFast(L,x),

    that inputs and outputs the SAME things as InterPol(L,x) but uses the ready-made Lagrange Interpolation Formula that you have just rediscovered. Compare the timing of

  5. Realize the drawback of equi-distance 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.

  6. Read about the so-called Runge phenomenon in section 3.2.1. (pp. 54-55) 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, Newton-Cotes Quadrature, and starting Gaussian Quadrature.

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

  1. Write a program
    ApplyNCg(f,x,n,a,b) that uses Newton-Cotes 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)
  2. 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.
  3. 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

  4. 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);
  5. 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:

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).

  1. [This problem is dedicated to my Patron Saint, Gian-Carlo Rota]
    The mapping
    P(x) -> int(P(x),x=0..1) ;
    is an example of LINEAR FUNCTIONAL (aka as Umbra) on the (infinite-dimensional) 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)

  2. Generalize procedure Pn(n,x) to write a procedure
    PnG(n,x,U,m)
    that inputs a non-neg. 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..n-1), 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.

  3. By using GuessRF, write a procedure
    GuessCoeff(n,i,U,m)
    that inputs a discrete variable n, a non-negative 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^(n-i) in PnG(n,x,U,m).
    What are

  4. [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 closed-form expression in terms of the symbols n and r, for the coefficient of x^(n-r) 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

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.1-5.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!

  1. [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*Pn(x)= Pn+1(x)+A(n)*Pn(x)+B(n)*Pn-1(x)
    Wrire an analogous procedure, BnM(M,n) to compute the first n terms of the sequence B(n) above.

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

  3. 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 n1-th term in the sequence of polynomials satisfying the recurrence
    x*Pn(x)= Pn+1(x)+A(n)*Pn(x)+B(n)*Pn-1(x)
    under the initial conditions
    P0(x)=1   P1(x)=x+a0  
    Hint: first rewrite the recurrence in a form expressing Pn(x) as a linear combination of Pn-1(x) and Pn-2(x) .

  4. 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: find (time permitting!, the last one may take too long), for each, the P1000(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

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!

  1. Do a second reading, and try to understand sections 5.1-5.3 of Herb Wilf's masterpiece.

  2. By plain brute force, or by the "method of bisection-searching", and using
    nops({op(NumerBifOrbF(f,x,g,x0,K1,K2,d1))})
    as an indicator of the probable orbit-size (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 bi-ai < eps, and such that when g=ai the generalized logistic map has period 2i-1 and when g=bi it has period 2i . For example,

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

    should yield something like   : [[2.99,3.01], [3.44,3.46]]   .
    What are

  3. [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*(1-x),x);
    should give 3.
    What are

  4. 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 power-series 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:

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!

  1. 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:

  2. Modify Zg(f,x,y,y0,N) to write a procedure
    PolAppx(f,x,y,y0,x1,N,eps):
    that inputs 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:

  3. [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!

  4. [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

    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:


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:

We also started numerical solutions of odes: C13a.txt , that contains procedures:

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!

  1. 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.

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

  3. 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 Runge-Kutta (RK4) gives you an approxination correct up to Di digits.

    What are

  4. Do a first reading of the Wikipedia entry on Runge-Kutta 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 Runge-Kutta 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.

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!

  1. 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.

  2. 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 weighted-average,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).

  3. Write a procedure, DelinZheng(lambda), that inputs a number lambda, and outputs the Runge-Kutta 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 so-called Butcher matrix. Make sure that DelinZheng(2)=[RK4].

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

  5. To prove that a given Runge-Kutta 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 semi-rigorously 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)()
  6. Write a procedure
    srProveB(B,d,M,K)
    that rigorously proves that B is NOT of order d if it returns false, and semi-rigorously 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 Runge-Kutta 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!

  1. Write a program FindRK(s,d) that inputs positive integers s and d and finds all Runge-Kutta methods with s steps and order d. For example,

    RK(4,4)

    should either give Delin-Zheng, or something even more general (i.e. with more than one parameter) for which Delin-Zheng 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 c-list 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 penny-pinching, and use the fact that the sum of the b's should be 1.

  2. [Challenge, I don't know the answer] Without peeking in the internet, find an order-5 Runge-Kutta method, using the present approach, possibly enhanced.
    [Added March 15, 2014, 4:38pm. According to L. Nick Trefethen's essay in p.604-615 of the "Princeton Companion of Mathematics" (W.T. Gowers, ed.), available from this secret url, there is no order-5 Runge-Kutta with 5 stages (i.e. s=5). Can you have Maple prove it rigorously? Can you find (w/o peeking) an order-5 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 Multi-step Adams-bashford methods; Implicit Runge-Kutte methods.

Contains the following procedures (and Butcher matrices RK4, GL)

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!

  1. 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 closed-form 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

    1. BestAdams(y,x,y,1)
    2. BestAdams(x+y,x,y,5)
    3. BestAdams(x+3*y,x,y,10)

  2. 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/4-1/6*sqrt(3)],[1/4+1/6*sqrt(3),1/4]],[1/2,1/2],[1/2-sqrt(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 warm-up.

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!

  1. [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(n-1)+L[2]*a(n-2)+ ... + L[-1]*a(n-nops(L))+f(n)

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

  3. 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 1-p, 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 i-1 dollars with probability 1-p.

  4. 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 i-th 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 i-1 dollars with probability 1-p, but you have wasted one round.

  5. Using the outputs of GR(1/2,N), and LGR(1/2,n), conjecture closed-form 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.

  6. Do first reading of chapter 3 of chapters 3-5 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 warm-up.

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!

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

  2. Using the above-mentioned 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 boundary-value recurrence equations.

  3. 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.

  4. 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)

  5. 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

    1. EmpiricalStability(x^4,x,0,1,[1/10,1/20,1/30,1/40]);
    2. EmpiricalStability(exp(x),x,0,1,[1/10,1/20,1/30,1/40]);
    3. EmpiricalStability(1/(3-x),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 1-Dimensional Heat Equation, following the wikipedia article of Finite Difference Methods Contains procedures:

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!

  1. Adapt procedure
    Heq2(U,x,h,k,x1,t1)
    to write procedure

    Heq3(U,x,h,k,x1,t1)

    that implements the Crank-Nicolson method descibed in the above-mentioned wikipedia article

  2. 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.

  3. [Challenge] By approximating the Laplace operator

    diff(u,x,x)+ diff(u,y,y)

    by the so-called five-point stencil

    (u(x+h,y)-2*u(x,y)+u(x-h,y))/h^2 + (u(x,y+k)-2*u(x,y)+u(x,y-k))/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 finite-difference approximation to u(x1,y1), using mesh size h in the x direction, and mesh-size k in the y-direction, 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

    1. Dirichlet(x^2+y^2,x,y,1/10,1/10,1/2,1/2)
    2. Dirichlet(x^2+y^2,x,y,1/20,1/20,1/2,1/2)
    3. 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