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:00noon-1:20pm)
• "Textbooks": Links to wikipedia and other internet sites
• Dr. Zeilberger's Office: Hill Center 704
• Dr. Zeilberger's E-mail: 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:30am-11:30am and by appointment

Description

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"

• (for novices only!) Frank Garvan's awesome Maple booklet ( part 1, part 2)
  Note: a few commands are no longer valid in today's Maple. The most important is that " has been replaced by %.
• L. Nick Trefethen's Definition of Numerical Analysis.
• L. Nick Trefethen's beautiful essay on Numerical Analysis, that appeared on pp. 604-615 of the "Princeton Companion of Mathematics" (W.T. Gowers, ed.), available from this secret url.
• the power method by anonymous.
• Wikipedia article about the power method.
• Herb Wilf's masterpiece, "Mathematics for the Physical Sciences"
• A very nicely written account on Chebyshev Expansions, that is chapter 3 of the book "Numerical Methods for Special Functions" by Amparo Gil, Javier Segura, and Nico Temme, published by SIAM in 2007.
• Dr. Z.'s Calc4 Handout 7
• The Wikipedia entry on Runge-Kutta methods
• wikipedia's article on Linear Multistep methods
• wikipedia article of Finite Difference Methods
• chapter 1 of István Faragó's great book

  [A more advanced (and detailed) treatment is given in this nice chapter, with unspecified author, on finite difference discretization of elliptic equations.]

• chapters 3-5 of Garrett Birkhoff's classic, "The Numerical Solution of Elliptic Equations".
• wikipedia's article on the Ritz method
• wikipedia's article on the Finite Element method
• Eric Weisstein's Wolfram Mathworld article on the Calculus of Variation
• Sections 11.5 and 12. 5 of the classic text on Numerical Analysis by R.L. Burden and J.D. Faires
• wikipedia entry on hearing the shape of a drum.
• Rodrigo Banuelos' wonderful lecture on the Weyl and Heat asymptotics.

Added March 23, 2014: Pick a project

Getting Started Homework (assigned Jan. 12, 2014):

How to submit your homework

• 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):

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)

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 d² and d^(log[2](3))= d^1.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) .

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(xn-alpha) 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)

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

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)

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

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)

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)

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!

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 one-dimensional world]

Homework for Thurs., Feb. 6 class (due Sunday Feb 9., 2014, 10:00pm)

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 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 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 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?

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:

• Hn(n): the determinant of the famous Hilbert n by n matrix (1/(i+j-1))(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+j-1+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 non-neg. 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)

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.

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:

• 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 ≤ n-1, 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 ≤ n-1, 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)

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

   • 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);

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

Programs done on Monday, Feb. 17, 2014

C8.txt , The Lagrange Interpolation formula, Newton-Cotes 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 (x-x1)*(x-x2)* ...*x(x-xn)

• 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 Newton-Cotes 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[n-1]*f((n-1)/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 Newton-Cotes 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 non-linear 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)

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
   • 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);

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?

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

• FindGaussSmart(n): Smart version of FindGauss(n) inputs a pos. integer n and finds, using non-linear 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*n-1 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)

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
   • 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 ?

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)

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ⁱ outputs the n-th 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_n-1(x)
  

Reading Homework for Monday Feb. 24 class (due 11:50am, Thurs. Feb. 27, 2014)

Homework for Monday Feb. 24 class (due Sunday March. 2, 2014, 10:00pm)

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*P_n(x)= P_n+1(x)+A(n)*P_n(x)+B(n)*P_n-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
   • 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)
     

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*P_n(x)= P_n+1(x)+A(n)*P_n(x)+B(n)*P_n-1(x)
   under the initial conditions
   P₀(x)=1   P₁(x)=x+a0  
   Hint: first rewrite the recurrence in a form expressing P_n(x) as a linear combination of P_n-1(x) and P_n-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:
   • 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₁₀₀₀(x). Convince yourself that it would take for ever to do it with PnM(M,1000,x) with the appropriate M.

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*(1-x) 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);
  

Homework for Thurs. Feb. 27 class (due Sunday March. 2, 2014, 10:00pm)

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 2^i-1 and when g=bi it has period 2ⁱ . 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

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

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

   • FirstBif(x^2*(1-x),x);
   • FirstBif(x*(1-x)^2,x);
   • FirstBif(x^2*(1-x)^2,x);
   • FirstBif(x^3*(1-x)^3,x);

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]

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:

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

Homework for Monday March 3, class (due Sunday March. 9, 2014, 10:00pm)

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:
   • 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];

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

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

   • (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 single-variable (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);

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 n-th partial sum of the first formula (for 128/Pi²) in Jesus Guillera's beautiful homepage.
• JesusPi(n): computes an approx. to Pi from the above appx. to 128/Pi^2

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

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

   • 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);

4. Do a first reading of the Wikipedia entry on Runge-Kutta methods.

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.

• RK4: NOT a procedure, but just an object, the Butcher representation of the classical Runge-Kutta
  [[],[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 Runge-Kutta with Butcher matrix a,b,c (see wiki)

• GRK(a,b,c,f,x,y,x0,y0,h,x1): The general Runge-Kutta 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 step-size h the output is pair of lists Xvalues, Yavlues

• [From C12.txt]:
  Z(f,x,y,y0,N): the N-th 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 Runge-Kutta with Butcher matrix a,b,c, BUT, only retaining powers of h up to R. Used for proving RIGOROUSLY, that the given Runge-Kutta 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 Runge-Kutta 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)

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

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.

• RP(x,y,M,d,N:=1): 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 N and M (note: N is optional and defaults to 1, but it's important to be able to test with negative coefficients sometimes)

• srProveB(B,d,M,K,N:=1): 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,N). (Note: N is optional and defaults to 1, but just in case you want to test negative coefficients, you can do so)

• DelinZheng(lambda): 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.

  [The above three procedures were written by Nathan Fox, taken from his completed homework 14]

• GButcher(a,b,s): a generic Butcher matrix of an explicit and coherent Runge-Kutta method with s stages output has the form [a,b,c] see above

• DiffRK(B,d,f,x,y,h): inputs a Butcher triplet [a,b,c] and a positive integer d, finds the difference between the actual beginning of the solution of the ODE y'(x)=f(x,y(x)), y(0)=0, and the output of the RK method of B up to order d in h

Homework for Thurs. March 13, class (due Sun. March. 16, 2014, 10:00pm)

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)

• 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) Gauss-Legendre Runge-Kutta method:
  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]:
• 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 Adams-Bashford method
    y(x+h) ~ y(x)+ h*add(L[i+1]*f(x-h*i,y[x-h*i]),i=0..nops(L)-1)
    (e.g. Euler is L=[1], 2-step Adams-B: 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 Adams-Bashdford d-step method (of order d)
• GIRK1(a,b,c,f,x,y,x0,y0,h): One step at the Implicit Runge-Kutte method with Butcher matrix a,b,c
• GIRK(a,b,c,f,x,y,x0,y0,h,x1): The general Runge-Kutta 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 step-size 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)

• wikipedia article of Finite Difference Methods
• chapter 1 of István Faragó's great book

Programming Homework for Monday, March 24, class (due Sun. March. 30, 2014, 10:00pm)

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.

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.

• GPi(N): the geometrical discrete Pi with mesh-size 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 mesh-size h
• Srec(L,Ini,n): The n-th term of the solution of the linear recurrence equation with constant coefficients
  a(n)=L[1]*a(n-1)+L[2]*a(n-2)+ ... + L[-1]*a(n-nops(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(n-1)+L[2]*a(n-2)+ ... + L[-1]*a(n-nops(L))
  a(0)=Ini[1], ..., a(nops(Ini)-1)=Ini[-1].
  a(N)=Fini[1], a(N-1)=Fini[2], ..
  

Homework for Thurs., March 27, class (due Sun. March. 30, 2014, 10:00pm)

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.

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.

• 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(n-1)+L[2]*a(n-2)+ ... + L[-1]*a(n-nops(L))
  a(0)=Ini[1], ..., a(nops(Ini)-1)=Ini[-1]
  a(N)=Fini[1], a(N-1)=Fini[2], ..
• P1ch(f,x,h,mu0,mu1): the exact solution of the 1-D 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 1-D 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(x-h)
  
• 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)

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]);

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:

• 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 1D-Heat 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): Floating-point 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)

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)

Programs done on Mon., April 7, 2014

• C20.txt: contains:
  • Dirichlet(f,x,y,h,k,x1,y1): appx. value at (x1,y1) of the sol. of the elliptic bvp

    u_{xx}+u_{yy}=f(x,y)

    with Dirichlet boundary conditions (i.e. 0 on the boundary) [code written by Cole Franks]

• C20a.txt: Looking for discrete versions of the largestpossible order for a differential operator with constant coefficients. Contains:
  • DisOp(Oper,D1,Sten,Sh,h) inputs a differential operator with constant coeff. phrased in terms of the differentiation operator D1, a set of integers, Sten, a symbol, Sh, for the shift operator

    Sh(u(x)):=u(x+h)

    a symbol h, for the mesh size

    Outputs: a recurrence operator (a Laurent poly in Sh) whose powers are in the Sten, with the highest possible order, along with that order

    For example

    DisOp(D1,D1,{0,1},Sh,h); should give

    (Sh-1)/h

  Homework for Monday, April 7, class (due Sun. April 13, 2014, 10:00pm)

  All homework should be uploaded to the secret directory that you gave me, as file hw20.txt. Please start the file with
  #hw20.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

     BVPexact(Oper,D1,f,x,a,b,Ini,Fini,x1)

     that uses dsolve to find the EXACT value of u(x1) where u(x) is a solution of the ODE Boundary value problem

     Oper(D1)u(x)=f(x) , f(a)=Ini[1], f'(a)=Ini[2], ...; f(b)=Fini[1], f'(b)=Fini[2], ...

     Oper is a differential operator written as a polynomial in D1, D1 is the symbol for the differentiation opeator, f is an expression in the variable x, that stands for the function f(x), a, b, are numbers with a What are [corrected April 9, 2014 (thanks to Andrew Lohr)]

     1. BVPexact(D1^2+1,D1,x^3,x,0,1,[1.5],[2.5],0.5);

     2. BVPexact(D1^4+3*D1^3+2*D1+1,D1,x^5+1,x,0,1,[1.0 , 2.0],[2.0 , 3.0 ],0.5);

  2. Do a first reading of section 11.5 that I handed out today, entitled "Rayleigh-Ritz method"
  
  Added April 14, 2014: See the students' nice Solutions to the 20th homework assignment

  Programs done on Thurs, April 10, 2014

  • C21.txt: contains:
    • BVP3(Oper,D1,f,x,h1,Ini,Fini,x1): inputs
      • Oper: a second-order const. coeff. diff. oper phrased in terms of the diff. oper. D1
      • an expression,f, in the cont. variable x
      • a symbol Sh for the shift operator u(x)->u(x+h) (h is a symbol)
      • a small NUMERIC mesh-size h1 (a reciprical of a positive integer)
      • x1 a number between 0 and 1 ( a multiple of h1)
      OUTPUTS:

      an appx. to u(x1) where u(x) is a solution BVP

      Oper(D1)u(x)=f(x), u(0)=Ini, u(1)=Fini
      using the Central 3-point discretation of Oper with mesh size h1

    • BVP5(Oper,D1,f,x,h1,Ini,Fini,x1): exactly as above, but using the the Central 5-point discretation of Oper

    • PDisOp1(Oper,D1,D2,Sten,Shx,Shy,h,ord): inputs a partial differential operator with constant coeff. a set of 2D lattice points, Sten, symbols for the shift operator Shx is: u(x,y)->u(x+h,y), and Shy is: u(x,y)->u(x,y+h), a symbol for h, the mesh size and a pos. intetger ord

      Outputs: a recurrence operator (a Laurent poly in Shx, Shy) whose powers are in the Sten For example
      PDisOp1(D1^2+D2^2,D1,D2,{[0,0],[-1,0],[1,0],[0,-1],[0,1]}, Shx,Shy, h,3) ;
      shoud give
      (1/Shx+1/Shy+Shx+Shy-4)/h^2

    • PDisOp(Oper,D1,D2,Sten,Shx,Shy,h,MaxOrd): same as above, but gives you the one with largest possible order.

      [Note: this is a slightly improved version of what we did in class]

    Homework for Thurs., April 10, class (due Sun. April 13, 2014, 10:00pm)

    All homework should be uploaded to the secret directory that you gave me, as file hw21.txt. Please start the file with
    #hw21.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

       BVPk(Oper,k,D1,f,x,h1,Ini,Fini,x1) that generalizes

       BVP3(Oper,D1,f,x,h1,Ini,Fini,x1) and BVP5(Oper,D1,f,x,h1,Ini,Fini,x1)

       by using, in the middle the stencil

       {-k,-(k-1), ..., k-1, k}

       and adjusting near the boundary such that
       BVPk(Oper,1,D1,f,x,h1,Ini,Fini,x1) and BVPk(Oper,2,D1,f,x,h1,Ini,Fini,x1)
       output the same thing as BVP3(Oper,D1,f,x,h1,Ini,Fini,x1) and BVP5(Oper,D1,f,x,h1,Ini,Fini,x1) respectively.
       

    2. Using
       BVPexact(Oper,D1,f,x,a,b,Ini,Fini,x1)
       that you wrote in hw20.txt, write a procedure
       Compare(Oper,D1,f,x,h1,Ini,Fini,x1)
       that outputs the quadruple
       [Exact Answer, Using k=1, Using k=2, Using k=3, Using k=4]
       with BVPk(Oper,k,D1,f,x,h1,Ini,Fini,x1) for k=1,2,3,4
       What are
       1. Compare(D1^2+4*D1+11,D1,x^5,x,1/20,1,2,1/2)
       2. Compare(D1^2-6*D1+7,D1,exp(x),x,1/20,2,3,1/2)

    3. [Challenge, but please try!]
       Using

       PDisOp(Oper,D1,D2,{seq(seq([i1,j1],i1=-2..2),j1=-2..2)},Shx,Shy,h,MaxOrd):

       as the Stencil for interior points, and appropriate analogs finite-difference versions for points near the boundary, adapt Cole Franks'
       Dirichlet(f,x,y,h,k,x1,y1):
       To write a procedure
       SuperDirichlet(f,x,y,h,k,x1,y1):

       Compare the outputs to an example where you can find the answer exactly (by starting with a function that vanishes on the boundary on the unit square (e.g. is a multiple of sin(Pi*x)*sin(Pi*y) by another function), and applying the Laplacian to it, and calling it f.

    4. Do a second reading of section 11.5 that I handed out Monday, entitled "Rayleigh-Ritz method"
    
    Added April 14, 2014: See the students' nice Solutions to the 21st homework assignment

    Programs done on Monday, April 14, 2014

    Today we started the Rayleigh-Ritz method, based on the Calculus of Variation following section 11.5 of Burden and Faires's book Numerical Analysis

    • C22.txt:
    contains:
    • EL(L,Y,Y1,y,t): implements the Euler-Lagrange Equation inputs an expression L in Y (stading for y(t)) Y1 (short for y'(t)) y and t are symbols

    • Eval1(Lx,Lf,x,x1): Given a partition of [0,1]=[x0,x1,x2,...,xn,x_{n+1}] and corresponding list of expressions Lf of n+1 expressions where the function between xi and x_{i+1} is Lf[i+1] EVALUATE it at x1

    • phii(Lx,x,i) : The function φ_i defined by Eq. (11.41)

    • TestF(c,Lx,x): a typical inear combination of φ_i

    Homework for Monday April 14, class (due Sun. April 20, 2014, 10:00pm)

    All homework should be uploaded to the secret directory that you gave me, as file hw22.txt. Please start the file with
    #hw22.txt, Date, and YOUR NAME

    If you do NOT wish your homework to be posted please write
    #Please do not post!

    1. 
       [Corrected Aprl 16, 2014, thanks to Frank Wagner]
       Write programs

       CubicSplineF(Lx,Lf,x), CubicSplineC(Lx,Lf,Lprimef, x)

       that inputs:

       • Lx and Lf, lists of numbers length 5 with Lx increasing,
       • Lprimef, a list of length 2 (only needed for CubicSplineC), that indicates the values of the derivate of the function at the endpoints x0 (alias Lx[1]), and x4 (alias Lx[5]),
       • x a symbol for the variable

       outputs a list of cubic polynomials in x of length 4 that produces a cubic spline according to the specs at the bottom of p. 595, for the free and clamped boundary conditions, respectively. The ouptut should be a list of cubic polynomials of length 4,

       [S0,S1,S2,S3]

       according to the notation there.

    2. Do a first reading of section 12.5 of Burden and Faires's book Numerical Analysis.
    
    Added April 21, 2014: See the students' nice Solutions to the 22nd homework assignment

    Programs done on Thurs, April 17, 2014

    • C23.txt:
    contains:
    • TestF(c,Lx,x): a typical linear combination of the phi_i, followed by the set of undetermined coefficients.

    • EvalFunct(p,q,f,x,Lx,Lu): inputs expressions p,q,f in the variable x, a partition, Lx, of [0,1], and a piece-wise function Lu, in x outputs the value of the functionl given in eq. (11.37) of the book, namely

      INT( p(x)*u'(x)^2+q(x)*u(x)^2-2*f(x)*u(x),x=0..1);

      RR(p,q,f,x,Lx): implements the Rayleigh-Ritz method to appx. the solution of the BVP ODE

      -(p(x)*y'(x))' + q(x)*y(x)=f(x), 0<=x<=1, u(0)=0, u(1)=0

    • [Under construction]
      BSpline(d,Lx,i), inputs a pos. integer d, and a list Lx=[x0, ..., x(d+1)] of numbers (increasing), outputs the B-spline of degree d on these points, i.e. it is 0 for xxd, and continuous to order d-1 (d+1)*(d+1) unknown numbers, (d+2)*d equations and f(Lx[i+1])=1

    Homework for Thurs. April 17, class (due Sun. April 20, 2014, 10:00pm)

    All homework should be uploaded to the secret directory that you gave me, as file hw23.txt. Please start the file with
    #hw23.txt, Date, and YOUR NAME

    If you do NOT wish your homework to be posted please write
    #Please do not post!

    1. Using Maple's command `piecewise', write a program

       Pw(Lf,x,Lx)

       that inputs a piecwise function in our format, where Lf is a list of expressions in the variable x, and Lx is a list of increasing numbers, and outputs it in Maple's notation. For example

       Pw([x,x^2],x,[0,1,2])

       should return

       piecewise(x<0,0, x>=0 and x<1,x, 1<=x and x<2, x^2, x>=2, 0)

    2. Write a procedure,

       Compare(p,q,v,x,x1)

       That inputs expressions p,q, v in x, and a number x1 between 0 and 1, then constructs u(x):=v(x)*x*(1-x) and

       f(x):= -(p(x)*u'(x))'+ q(x)*u(x)

       (so we know a priori the exact solution of the BVP

       -(p(x)*y'(x))'+ q(x)*y(x)=f(x) , 0<=x<=1, y(0)=0, y(1)=0

       namely u(x)

       and then uses

       RR(p,q,f,x,[seq(i/N,i=0..N)]):

       with N=10,20,30,40, and evlautes them at x1. The output is a list of length 5 with the last entry being the exact value.

       What are

       1. Compare(x^2,1+x^3,sin(x),x,.513);
       2. Compare(x^2,1+x^3,x^7+1,x,.513);
       3. Compare(x^2+3,sin(x),exp(x),x,.542);

    3. [Corrected April 20, 2014, 3:32pm, thanks to Anthony Zaleski]

       Finish up procedure

       BSpline(d,Lx,i)

       Verify that

       BSpline(3,[-2,-1,0,1,2],2);

       agrees with the piecewise function S(x) given by Eq. (11.43) on page 596 of Burden and Faires's book Numerical Analysis.

    4. Write an analog of RR(p,q,f,x,Lx), calling it

       RR3(p,q,f,x,Lx)

       where the phi_i(x) are replaced by those defined on p. 597, using the cubic B-spline, S(x)

    5. Do a second reading of section 12.5 of Burden and Faires's book Numerical Analysis.
    
    Added April 21, 2014: See the students' nice Solutions to the 23rd homework assignment

    Programs done on Monday, April 21, 2014

    • C24.txt:
    contains:
    • ET1(T,A,x,y): inputs a triangle T (given as a set of three vertices), and A, a member of T, and variables x,y outputs the (uniqe) linear expression a+b*x+c*y that is 1 at A and 0 at the other two vertices

    • Le(A,B,x,y): inputs two points A and B and outputs the line through them

    • InsT(T,x,y): the inside of the triangle T (modified from class by relacing the inequalities by strict ones)

    • ET(T,A,x,y): inputs a triangle T (given as a set of three vertices), and A, a member of T, and variables x,y outputs the (uniqe) piecewise expression that equals something like a+b*x+c*y INSIDE the triangle, such that is 1 at A and 0 at the other two vertices and 0 outside the triangle

    • Tg(N): tirangulates into right-angled triangles of edge-side 1 the square [0,N]x[0,N]

    • TestF(T,c,x,y): the generic linear combination of elements that defines a piece-wise continuous function on the domain that is the union of the insides of the triangles of the set T

    Homework for Monday April 21, class (due Sun. April 27, 2014, 10:00pm)

    All homework should be uploaded to the secret directory that you gave me, as file hw24.txt. Please start the file with
    #hw24.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, InsP(P,x,y), where P is a list of vertices, describing, in counterclockwise way, the vertices of a polygon assumed to be convex, and outputs the conditions for a point to be inside.

    2. Write a program TgT(N), that inputs a positive integer N and outputs a natural triangulation of the equilateral triangle whose vertices are at
       [0,0],[1,0],[1/2,sqrt(3)/2]
       into equilateral triangles of side-length 1/N. What are
       1. TgT(2)
       2. TgT(4)

    3. Write a program TgC(N), that inputs a positive integer N and outputs a natural triangulation APPROXIMATING the unit circle, centered at the origin, into right-angled triangles of side-length 1/N What are
       1. TgC(2)
       2. TgC(8)
    
    Added April 28, 2014: See the students' nice Solutions to the 24th homework assignment

    Programs done on Thurs., April 24, 2014

    • C25.txt, continuing the Finite Element Method.

    Contains procedures:

    • TestF(T,c,x,y): the generic linear combination of elements that defines a piece-wise continuous function on the domain that is the union of the insides of the triangles of the set T, expressed in OUR data structure as a list of pairs [triangle, expression in x,y]

    • IntTr(f,x,y,t): double-integral of f over the interior of the triangle t [discovered by Nathaniel Shar]

    • Iw(F,p,q,r,f,x,y): inputs a piece-wise function F (given a list of pairs [t,ExpIn(x,y)] and expressions p,q,r,f in the variables x,y, and a triangle t, finds I[w] as given in (12.44) of the book

    • FEM(p,q,r,f,g,x,y,N): the finite-element apprx. solution to the elliptic PDE Dirichlet BVP(12.41), in the region

      { (x,y) | 0 ≤ x ≤ N , 0 ≤ y ≤ N }

      with u=g on the boundary using the triangulation Tg(N)

    • [Finished after class]
      EvalF1(F,x,y,pt): evaluate the point pt in the piecewise function F

    Homework for Thurs. April 24, class (due Sun. April 27, 2014, 10:00pm)

    All homework should be uploaded to the secret directory that you gave me, as file hw25.txt. Please start the file with
    #hw25.txt, Date, and YOUR NAME

    If you do NOT wish your homework to be posted please write
    #Please do not post!

    1. Work hard on your chosen project. The first draft is due May 6, 2014, 10:00pm. It is understood that many of you would continue to develop it in the summer and beyond, possibly turning it into a publishable article.

    2. Modify procedure Tg(N) to write a procedure
       Tgh(N,h)
       where h is a reciprocal of a positive integer, and outputs an analogous triangulation with right-angled triangles of base and height h. Note that Tgh(N,1) and Tg(N)should output the same thing for every N.

    3. Modify procedure
       FEM(p,q,r,f,g,x,y,N)
       done in class, to write procedure
       FEMh(p,q,r,f,g,x,y,N,h)
       that uses the triangulation Tgh(N,h) that you wrote above, don't forget to redefine the members of the boundary, B.

    4. Wtite a procedure
       CheckFEM(u,p,q,r,x,y,pt)
       that starts with a solution, u (an expression of x,y), expressions p,q,r as in (12.41) of the book, computes f, and then computes
       FEMh(p,q,r,f,u,x,y,1,h)
       (note that g is u, since we already know the solution and it better be equal to g on the boundary) at
       h=1/4,1/8,1/12,1/16
       and plug-it in at the point p. The output should be a list of four numbers, followed by the EXACT value of u at (x,y)=pt.

       What are

       1. CheckFEM(x+3*y+10,x+y+3,-3-x^2*y,x^2+y,x,y,[0.32,0.52]);

       2. CheckFEM(x^2*y+7*x*y+2,x^2+y+5,-4*x-y,x^3+3*y,x,y,[0.32,0.52]);

    
    Added April 28, 2014: See the students' nice Solutions to the 25th homework assignment

    Programs done on Monday, April 28, 2014

    • C26.txt, Gambler's ruin in one dimension

    Contains procedures:

    • Die1: inputs p rational between 0 and 1 and outputs 1 with prob. p and 0 with prob. 1-p

    • GR(i,N,p): simulates a gambler's ruin session staring with i dollars, ending with 0 or N dollars and with prob. of winning a dollar =p

      ExpD(N,p): Let f(i) be the expected duration of the game starting at i, outputs the vector [f(1),f(2),..., f(N-1)] of expected durations with prob. p of winning one dollar
      f(i)=p*f(i+1)+(1-p)*f(i-1)+1, i=1..N-1, f(0)=0, f(N)=0

    • PGF(N,p,t): Let f(i,t) be the prob. generating function of the a rational function in t whose Maclaurin expansion's coeff. of t^j gives the exact prob. of finishing in j steps game starting at i, outputs the vector [f(1),f(2),..., f(N-1)] of expected durations with prob. p of winning one dollar f(i,t)=t*(p*f(i+1,t)+(1-p)*f(i-1,t)), i=1..N-1, f(0)=0, f(N)=0

    • PGFf(N,p,t): fast version; of PGF(N,p,t)

    • [Added after class]
      GuessP(L,x): guesses a polynomial expression in x for the list L such that L[i]=P(i)

    Homework for Monday April 28, class (due Sun. May 4, 2014, 10:00pm)

    All homework should be uploaded to the secret directory that you gave me, as file hw26.txt. Please start the file with
    #hw26.txt, Date, and YOUR NAME

    If you do NOT wish your homework to be posted please write
    #Please do not post!

    1. Work hard on your chosen project. The first draft is due May 6, 2014, 10:00pm. It is understood that many of you would continue to develop it in the summer and beyond, possibly turning it into a publishable article.

    2. Write a two-dimensional analog of ExpD(N,1/2), call it ExpD2(M,N), that inputs positive integers, and outputs a list of lists, let's call it L, such that

       L[i][j] (1 ≤ i ≤ M-1   ,   1 ≤ j ≤ N-1 )

       is the expected number of unit steps left for a random walker who starts at [i,j] and moves, in a unit step, with equal probability (1/4) to one of its immediate neighbors

       [i-1,j],[i+1,j],[i,j-1],[i,j+1]

       and ends the walk when reaching one of the edges i=0, i=M, j=0, j=N

    
    Added May 5, 2014: See the students' nice Solutions to the 26th homework assignment

    Programs done on Thurs., May 1, 2014

    • C27.txt, More Gambler's ruin in one dimension, giving explicit expressions for the first few moments about the mean of the r.v. "duration of the game".

    Contains (in addition to those of C26.txt that are also included), procedures

    • PGFam(N,p,t): the list of probability generating functions in the variable t, for the r.v. "duration of a gambler's ruin game with exit capital N" about the mean, with starting capital i (i from 1 to N-1)

    • Moms(N,p,R): the list of length N-1 whose ith item is the R-th moment of the above r.v.

    • Momsi(N,R,i): The explicit polynomial expression in the symbol i for the R-th moment of the r.v. remainign duration of the FAIR game with numeric largest capital N.

    • MomsiN(N,R,i): a guessed polynomial in symbolic i and N for the R-th moment of the r.v. "duration of a game staring at i dollars and ending at either 0 or N, with a FAIR coin

    • AStmomsiN(R,x): The leading term of the asympotic standarized R-th moment if you start with x*N , the output is an expression in x.

    Homework for Thurs. May 1, class (due Sun. May 4, 2014, 10:00pm)

    All homework should be uploaded to the secret directory that you gave me, as file hw27.txt. Please start the file with
    #hw27.txt, Date, and YOUR NAME

    If you do NOT wish your homework to be posted please write
    #Please do not post!

    1. Work hard on your chosen project. The first draft is due May 6, 2014, 10:00pm. It is understood that many of you would continue to develop it in the summer and beyond, possibly turning it into a publishable article.

    2. Modify procedure

       ExpD(N,p)

       to write a procedure

       ProbW(N,p)

       that inputs a positive integer N, a number p between 0 and 1, and outputs a list of length N-1 whose i-th entry tells you the probability of exiting the casino a winner, if you start out with i dollars, and must leave as soon as you are either broke or have N dollars.

       Use this procedure to find the probability of doubling your fortune if you currently have 50 dollars, rather than losing it with

       1. p=1/2
       2. p=17/35
       3. p=49/100
       convince yourself that it is a bad idea to gamble!

    3. If you have not done it already, fill-out the on-line teaching evaluation form, and don't even think of entering anything below 5.
    
    Added May 5, 2014: See the students' nice Solutions to the 27th homework assignment

    Lecture of May 5, 2014

    Today, the last class of this exciting course, was the most exciting. We had the honor and pleasure of a lecture by the one and only Neil Sloane, who talked, amomg other things on this sequences

    Homework for Monday May 5, class (due Tue. May 6, 2014, 10:00pm)

    1. Work hard on your chosen project. The first draft is due May 6, 2014, 10:00pm. It is understood that many of you would continue to develop it in the summer and beyond, possibly turning it into a publishable article.

    2. Have a great summer!

    Added May 7, 2014: Look at the students' great projects.

    ---

    Added May 19, 2014: Here are the students' evaluations.

    ---

    Dr. Z.'s teaching page