18.330 Introduction to Numerical Analysis
Spring 2015
Problem Set 10
Due: May 12, 2015, in class
Problem 1: Arbitrary-Precision Arithmetic. Use FFT techniques to implement computer subroutines for arbitrary-precision integer addition and multiplication. (Actually, addition is easy; only multiplication requires FFT techniques.1 ) Use your subroutine to compute the 100th Catalan number C100 ,
where
C0 = 1,
C1 = 1,
C2 = 2,
C3 = 5,
Cn+1 =
n
X
Ck Cn−k .
k=0
Note that C100 has 57 digits. For a list of 66 different ways to interpret the Catalan numbers, see Richard Stanley’s book Enumerative Combinatorics, Volume
2, Exercise 6.19. (And you thought my problem sets were long.)
Problem 2: Numerical differentiation by FFT.
In previous problem sets you investigated methods of numerical differentiation—that is, methods of estimating the derivative of a function f (x) given
only discrete samples of the function at discrete values of x. In this problem
you will explore the use of the FFT as a technique for numerical differentiation.
In this problem, let fn ≡ f (nh), n = 0, · · · , N − 1 be samples of f (x) at N
evenly-spaced data points.
(a) Write a computer subroutine that accepts a vector of function samples {fn }
and a number h and uses the (second-order) centered-difference stencil to
1 The easiest way to do this is to represent the decimal digits of an arbitrary-size integer
as an array of ordinary integers, with each array element lying between 0 and 9. Adding two
numbers then amounts to a simple elementwise addition of the digit arrays, with the slight
complication that the result of this procedure may involve “digits” that are greater than 10;
remedying this simply involves a “rectification” step (equivalent to “carrying” in elementaryschool arithmetic) in which each digit is reduced to a number between 0 and 9, possibly
by augmenting a neighboring digit. Multiplication corresponds to a discrete convolution of
the digit arrays (which you will here implement using FFT techniques) followed again by a
rectification step. Don’t forget zero padding!
1
compute and return a vector of derivative values {fn0 }. For differentiating
at the endpoints, your routine may assume that the underlying function
is periodic, so that f0 = fN and fN +1 = f1 . (Note: this subroutine should
be approximately one line long.)
(b) Now write a subroutine that differentiates using FFT techniques. This
subroutine will again accept a vector of function samples and a spacing
h. It will compute the DFT of the function samples and use the results
to construct a trigonometric interpolant passing through the data points.
[In testing your routine, you may find it useful to plot the trigonometric
interpolant you construct and confirm that it does run through all the data
points.] Then your program will differentiate this interpolant to compute
the values fn0 that it will return in the output vector.
(c) Compare the accuracy of your subroutines in parts (a) and (b) by testing
them on the function f (x) = esin(x) . That is, construct a vector of 25
evenly-spaced samples of this function at points x between 0 and 2π.
Also evaluate the derivative f 0 (x) analytically and use this expression to
construct a vector of exact derivative values at the same evenly-spaced
x points. Then plot, as a function of x, the error in the finite-difference
derivative and the error in the FFT derivative as computed using your
solutions to parts (b) and (c).
(d) For the function f (x) considered above, how would you expect the (i) computational cost (CPU time), (ii) accuracy of the FFT-based differentiation
to vary with the number of sample points N (where N = 25 in the above
example)? How do these predictions compare with the cost and accuracy
scaling of the finite-difference approach of part (a)? For the accuracy in
particular, explain why the FFT algorithm exhibits the behavior it does.
(e) Now consider using the FFT-based algorithm to differentiate the function
g(x) = ex over the same interval. How would you expect the (i) cost, (ii)
accuracy of the FFT-based differentiation algorithm to vary with N in
this case? [Just state your predictions; you don’t need to implement the
calculation in this case, although you may find it instructive to do so.]
Problem 3: Chebyshev differentiation and boundary-value problems.
In this problem you will revisit the boundary-value problem you investigated
on PSet 4. Recall that this problem was
d2 x
+ ω02 x = cos(t2 ),
dt2
x(0) = 0.9,
x(10) = −2.1,
ω0 = 3.
(1)
In PSet 4 you solved this problem using finite-difference techniques, which exhibit algebraic convergence. Here you will study the same problem using a
spectral method and compare the rates of convergence.
2
(a) Let f be a vector of length N + 1 whose entries are samples of a function f
at the Chebyshev points, i.e.
nπ fn = f (tn ),
tn = cos
N
for n = 0, 1, · · · , N . Similarly, let f 0 be the vector whose entries are
samples of the derivative of f :
nπ fn0 = f 0 (tn ),
tn = cos
.
N
Write a computer program that accepts an integer argument N and computes the (N + 1) × (N + 1) matrix DN that operates on f to yield f 0 .
(b) Write a computer program that uses your answer to Part (a) to solve the
BVP (1). [To ensure that your code is correct, you may wish to plot the
solution x(t) that it predicts and compare against the solution obtained in
Problem 2 of PSet 4.] Now compare the values computed for the particular
value x(5) by your Chebyshev code and by your finite-difference code.
Plot, as a function of the dimension of the linear system solved in each
case, the error in this quantity as computed by the two solution methods.
3
Extra credit (15%). In linguistics, a palindrome is a word or phrase that
remains invariant under left-right reversal. Examples include
• A Santa dog lived as a devil god at NASA.
• No, Mel Gibson is a casino’s big lemon.
• Yo, banana boy!
On the other hand, in mathematics a palindrome is an integer that remains
invariant under a left-to-right reversal of its decimal digits. Examples include
the integers 2552 and 42677624.
Now consider the following algorithm for producing a palindromic number
given an arbitrary initial integer n:
1. Set n = n + REVERSE(n).
Here REVERSE(n) is the integer whose decimal digits are the left-to-right
reversal of the decimal digits of n: for example, REVERSE(259) = 952.
2. Repeat step 1 until n is a palindrome, i.e. we have n=REVERSE(n).
For example, given the initial value n=194, the algorithm terminates in 3 iterations:
194 + 491 = 685
685 + 586 = 1271
1271 + 1721 = 2992.
A longstanding conjecture holds that this algorithm terminates for any starting
value, and evidence supporting the conjecture comes from the fact that, indeed,
for most starting values the iteration converges in a small number (∼ 20) of
iterations. However, confidence in the conjecture is shaken by the existence of
rare starting values for which the iteration seems to take an exceedingly long
time to converge. For example, there are 249 integers less than 10,000 for which
the iteration fails to converge within 100 iterations.
Using your implementation of arbitrary-precision arithmetic, plot the number of iterations required for the above iteration to converge—and the number
of digits in the final palindromic number—for starting values of n in the range
10 < n < 300. Do you observe any regularity in the pattern here? Identify at
least two integers for which the iteration fails to converge within 100 iterations.
References. A lovely discussion on mathematical palindromes may be
found in the article “Backward run numbers, letters, words, and sentences until
boggles the mind,” by Martin Gardner in the August 1970 issue of Scientific
American, http://www.nature.com/scientificamerican/journal/v223/n2/
pdf/scientificamerican0870-110.pdf.
Meanwhile, a classic meditation on the crucial psychosomatic benefits of
linguistic palindromes (which includes, among other gems, the instant classic
Lana C. Ladaug was I ere I saw Guadalcanal ) is “Ainmosni,” by Roger Angell
in the May 31, 1969 issue of the New Yorker : http://archives.newyorker.
com/?i=1969-05-31#folio=032.
4