Exercises in Numerical Analysis

Computing Errors

1. Let T be the transformation of the interval [0, 1) defined in class.

i. Which points x ∈ [0, 1) do have the property that T^{n}(x) = 0 for all sufficiently large n?

ii. Prove that, if x is not such a point, then we cannot even possibly have

T^{n}(x) ⟶_{n→∞} 0.

2. (requires Probability Theory) Use Borel-Cantelli Lemma to prove that, for l = 0, 1, … 2^{k} − 1, if x is chosen randomly (uniformly) in the interval [0, 1), then there is probability 1 that the interval [l/2^{k}, (l+1)/2^{k}) will contain infinitely many terms of the sequence (T^{n}(x))_{n=0}^{∞}. Conclude that, moreover, if x is chosen randomly in [0, 1), then there is probability 1 that every interval [a, b) ⊆ [0, 1) of positive length will contain infinitely many terms of that sequence.

3. i. Use the infinite binary expansion of real numbers to find a point x for which the sequence T^{n}(x) is asymptotically 30% of the time in the subinterval [0, 1/2) (and 70% of the time in the subinterval [1/2, 1)).

ii. Does there exist a point x such that T^{n}(x) lies in the subinterval [0, 1/4) asymptotically 40% of the time, in [1/4, 1/2) – 30% of the time, in [1/2, 3/4) – 20% of the time and in [3/4, 1) – 10% of the time?

4. Consider the family of transformations T_{a}, a = 2, 3, …, of the interval [0, 1), defined by T_{a}(x) = {ax}, where {t} denotes the fractional part of a real number t.

i. Given a transformation S of an arbitrary space X (namely, S is a function from X to itself), a point x ∈ X is eventually periodic under S if there exist positive integers m < n such that S^{n}(x) = S^{m}(x) (or, equivalently, if the orbit of x under S is finite). Which points of [0, 1) are eventually periodic under T_{a}?

ii. In the setup of the preceding part, a point x ∈ X is periodic under S if S^{n}(x) = x for some n ∈ N. Perform a few experiments, and conjecture which points of [0, 1) are periodic under T_{a}.

5. Let T_{a} be as in the preceding question.

i. Test how T_{a}-orbits look like on the computer for odd numbers a. Do you think the computer results are closer to the truth than they are for a = 2?

ii. Now consider even values of a. You can easily verify that, for some of them, the orbits received on the computer contain much less non-zero terms than in the case of a = 2, whereas for others these orbits contain just a few less such terms. Explain this dependence on a.

6. Given any transformation T of an interval (or a more general space), one may consider orbits {T^{n}(x) : n = 0, 1, 2, …} of points x. In the case considered in class, we saw that, if the transformation is defined as multiplication by 2 modulo 1, then even a very small error in the initial determination of x causes the computed values of T^{n}(x) to be quite far from their true values already for medium-sized n. Can you suggest a sufficient condition on T under which the initial error causes only small errors in the computation of T^{n}(x) even for large values of n? (Hint: Formulate your condition in terms of the relation between the size of |T (x_{2}) − T (x_{1})| and that of |x_{2} − x_{1}| for various pairs of points (x_{1}, x_{2}).)

Computer Representation of Numbers

7. i. Write a MATLAB function which, given a vector containing the digits in the decimal expansion of a non-negative integer L, provides the binary expansion of L.

ii. Write a MATLAB function which, given a non-negative integer, provides its binary expansion.

8. Consider the following MATLAB function, whose primary purpose is to calculate (the decimal representation of) a number, given its digits in some other base:

function n = convert(digits, base)

initial flops = flops;

n = 0;

for i = 1 : length(digits) n = base ∗ n + digits(i);

end

additional flops = flops − initial flops end; Explain why the function works properly, and what are its advantages with

respect to the function one would obtain by implementing in a straightforward manner the algorithm hinted in class.

9. i. Write a MATLAB function which, given an integer b ≥ 2 and two non-negative integers m and n in base b, calculates m+n. Assume that m and n are given as strings (of arbitrary length) of digits in base b, where, for simplicity, if b > 10 then the digits representing 10, 11, … are the characters right after ′9′ in ASCII. (That is, 10 is represented by ’:’, 11 – by ’;’, and so forth.)

ii. Do the same for multiplication instead of addition.

10. It is well known (prove!) that a positive integer L is divisible by 11 if and only if the alternating sum of the digits in its decimal expansion is divisible by 11. Namely, if the digits are a_{0}, a_{1}, … a_{k}, then the number is divisible by 11 if and only if a_{0} − a_{1} + a_{2} − … + (−1)^{k}a_{k} is divisible by 11. For which number, instead of 11, do we have the analogous criterion for divisibility when using binary expansion?

11. In base 10, one may decide whether or not a given number is divisible by 2^{k} or by 5^{k} by considering its k last digits. More generally, divisibility by 2^{k}5^{l} depends on the last max(k, l) digits only. For which numbers a can you decide divisibility by a, for a given number represented in base 2, by considering the k last digits of the number only?

12. It is quite cumbersome to compute the binary representation of numbers given in decimal, and vice versa. For which pairs of bases d_{1} and d_{2} is it easy to convert from one to the other?

13. Find the fixed point representation and the floating point representation of the following numbers: 511.75, −19.3125, \frac{3}{5}, −\frac{1}{7}.

14. Write down all positive floating point numbers for a system using 8 bits for representing numbers, out of which 3 bits are for the exponent field and 4 for the fraction field.

15. i. Calculate the sum of all positive numbers in a fixed point system using i bits for the integral part field and f for the fractional part field.

ii. Do the same for floating point numbers, using e bits for the exponent field and f for the fraction field.

16. Let η denote the minimal floating point number. Find the ratio \frac{η}{ε}, assuming that the system uses e bits for the exponent field and f for the fraction field.

17. i. Find the least positive integer which is not a floating point number. ii. Find the least positive integer a such that both a and a + 1 are not

floating point numbers.

iii. Same if all three numbers a, a + 1 and a + 2 are required not to be floating point numbers.

iv. Same for four consecutive numbers.

18. (riddle) Find the next term in the following finite sequence:

10, 11, 12, 13, 14, 15, 16, 17, 20, 22, 24, 31, 100, 121, ? .

The IEEE Standard; Rounding

19. Estimate the product of all positive normalized floating point numbers.

20. Performing appropriate calculations with MATLAB, show that the system:

i. employs the IEEE standard with double precision;

ii. has round-to-nearest as its rounding mode; iii. performs correctly rounded arithmetic.

21. Denote by max the largest floating point number in double precision. Use MATLAB to find the maximal number y which, when added to max, does not give the answer Inf. Explain!

22. i. For what range of small positive numbers η will the result of the expression 1 − cos(η) be 0 in double precision? Why?

ii. Same for e^{η^{3}} − 1.

23. i. Verify in MATLAB that there is a distinction between 0 and −0.

ii. Which arithmetic operations have −0 as their true result (even if they are presented as 0)?

24. Suppose two n × n matrices A and B are chosen with each entry of the two matrices being either 0 or −0, with equal probabilities. Let C = AB. We want to find the average number of 0 entries in C and the probability that C consists only of 0 entries.

i. Write a MATLAB function which calculates these two quantities exactly by going over all possible matrices A and B.

ii. Write a MATLAB function which estimates the two quantities by selecting randomly some required number of matrices A and B.

iii. (requires Probability Theory) Calculate exactly the first quantity.

25. Is it possible that, for two normalized positive numbers a and b, both of the quantities

c = \frac{a}{\sqrt{a^{2} + b^{2}}} , d = \frac{b}{\sqrt{a^{2} + b^{2}}}

(which are the cosine and sine of the argument of the complex number a + bi, respectively) will come out as Inf?

26. i. Find the largest positive integer n for which (3/2)^{n} is a floating point number in double precision.

ii. Write a MATLAB function which, given a positive integer n, provides the exact binary expansion of (3/2)^{n}. (The function should work, at least in principle, for any n.)

27. Same as the preceding question for n! instead of (3/2)^{n}.

28. Is the operation ⊕ commutative? associative?

29. Consider the example given in class of orbits of the transformation T_{2}. i. How many steps at most are needed for a floating point number in

[0, 1) to reach 0 under successive applications of T_{2}?

ii. (requires Probability Theory) Estimate the average number of iterations needed until a randomly selected number reaches 0.

Nonlinear Equations

30. i. Find a zero of the function e^{x} − 2 cos x in the interval [0, 1] with an error not exceeding 0.01 by means of the Bisection method, using only a calculator. Does the function have additional zeros?

ii. Same for the function e^{x} sin x − 1.

31. It is required to find all the zeros of the function f defined by

f(x) = xe^{2x} − 10xe^{x} − 5e^{x} + 50 .

i. Show that f has a single zero in the interval [1, 2], another single zero in the interval [2, 3], and no additional zeros.

ii. How many flops are required to find each of the zeros with an error not exceeding 10^{−6} if we use the Bisection method in MATLAB?

iii. Simplify the form of f to obtain a better result (still using the Bisection method when required).

32. Same as the preceding question for the function

f(x) = xe^{2x} sin x − 2xe^{x} − 5e^{x} sin 2x + 20 cos x ,

where the intervals [1, 2] and [2, 3] are replaced by [0, 1] and [1, 2], respectively.

33. i. Suppose instead of the Bisection method we use “trisection”. Namely, at each step we divide the current interval into three equal subintervals, and check which of them contains a zero of the given function. Assuming that the function contains a single zero in the interval [a, b], how fast would you expect the method to converge “on the average”? Compare the result with that of the Bisection method. (The analysis should be rather heuristic. Assume that, at each step, the zero of the function may belong to each of the three subintervals with equal likelyhoods.)

ii. Same as the first part for the following variant of the Bisection method: At each step the current interval [a′, b′] is bisected into two subintervals [a′, (1 − α)a′ + αb′] and [(1 − α)a′ + αb′, b′], where 0 < α < 1 is fixed. How does your answer vary with α?

34. Suggest a few ways of rewriting Problem 24.i as one of locating fixed points of an appropriately chosen function.

35. i. Let f be twice differentiable in the interval [a, b]. Assume f has at least one zero at the interval. Show that Newton’s Method leads within one step to a zero of f for every starting point x_{0} ∈ [a, b] if and only if f is linear.

ii. Construct a function having an infinite number of zeros (ξ_{n})_{n=1}^{∞}, so that each ξ_{n} has the following property: There exist infinitely many points (θ_{n,i})_{i=1}^{∞} such that Newton’s Method, if started at any θ_{n,i}, leads within one step to ξ_{n}.

36. i. Find a family of cubic polynomials for which Newton’s Method, if started at an appropriate point, yields a periodic orbit.

ii. (requires Differential Equations) Find a function f defined over the whole real line with the following properties: 1) f has a single zero. 2) f is differentiable everywhere, except perhaps at its zero. 3) Newton’s Method, with any initial point except for the zero of f, yields a periodic orbit.

37. Consider the following functions from R to itself. Characterize the fixed points of each, and determine which points are attracted to each of these under successive applications of the function. If possible, try to find the rate of convergence of those points attracted to the fixed points.

a. g(x) = \frac{2}{3}x + 7 . b. g(x) = 2x + 5 . c. g(x) = x + x^{3} . d. g(x) = x − x^{3} . e. g(x) = x + x^{2} . f. g(x) = sin x . g. g(x) = cos x .

h. g(x) = \{ \begin{matrix}x + x^{3} sin \frac{1}{x}, & x ≠ 0 , \\ 0, & x = 0 . \end{matrix}

38. The same as the preceding question for the following functions from R^{2} to itself: a. g(x, y) = (x/2, y/3) . b. g(x, y) = (x/2 + 1, y/3 + 2) . c. g(x, y) = (2x, 3y) . d. g(x, y) = (x/2 + 1, 3y − 2) . e. g(x, y) = (x/2 + y, y/2) . f. g(x, y) = (x + y, −2x + 4y) . g. g(x, y) = (x cos α − y sin α, x sin α + y cos α), (α ∈ R) . h. g(x, y) = (x^{3}y^{4}, x^{−4}y^{3}) (here take R_{+}^{2} as your domain.)

39. For the following functions all zeros can be explicitly found. Find the biggest interval you can around each zero such that, starting Newton’s Method at a point belonging to this interval, we obtain a sequence converging to that zero. Determine the rate of convergence in each case.

a. f(x) = e^{x} − a . b. f(x) = x^{2} − 8x + 15 . c. f(x) = (x − a)(x − b) . d. f(x) = (x − 1)(x − 2)(x − 3) . e. f(x) = (x − 1)^{k} . f. f(x) = sin x . g. f(x) = e^{x} − (1 + x + x^{2}/2) . h. f(x) = e^{x} − 1 − sin x , x ∈ [−\frac{π}{2} , \frac{π}{2} ] . i. f(x) = tgx .

40. Show that, if ξ is a root of order k of f, then, defining g by

g(x) = x − k \frac{f(x)}{f′(x)} ,

we have g′(ξ) = 0.

41. i. Suppose Newton’s Method is applied to the function f(x) = x^{k}. Show that, the larger is k, the slower is the convergence.

ii. Consider the function defined by:

f(x) = \{ \begin{matrix}e^{−1/x^{2}}, & x ≠ 0, \\ 0, & x = 0 . \end{matrix}

Find the rate of convergence of Newton’s Method for this function. Explain.

42. How does Newton’s Method work for the function f defined by:

f(x) = \{ \begin{matrix}\sqrt{x}, & x ≥ 0, \\ −\sqrt{−x}, & x < 0 . \end{matrix}

Explain how this corresponds to the theory.

43. Consider the following two iterative formulas, designed to solve the equation x^{2} = a:

\begin{matrix}x_{n+1} = & \frac{x_{n} + a/x_{n}}{2} ,  n = 0, 1, 2, …, \\ x_{n+1} = & \frac{x_{n} + a/x_{n}}{2} − \frac{(x_{n} − a/x_{n})^{2}}{8x_{n}} ,  n = 0, 1, 2, …. \end{matrix}

i. Show that the first method converges quadratically. ii. Show that the second method is a third-order method. iii. Write a MATLAB function implementing the two methods, and

compare their performance.

44. Consider the function defined by:

f(x) = \{ \begin{matrix}x^{4} sin \frac{1}{x}, & x ≠ 0 \\ 0, & x = 0 . \end{matrix}

Does there exist an interval around 0 such that, starting Newton’s Method at any point x_{0} of this interval with f(x_{0}) ≠ 0 and f′(x_{0}) ≠ 0, we obtain an orbit converging to 0?

45. i. Suppose f is defined on an interval [a, b] and has the property that, starting the Secant Method at any two points within the interval, we are led to a zero of f within a single step. Prove that f is necessarily linear.

ii. Construct an infinitely differentiable non-linear function f having the following property: For every point x_{0} there exist infinitely many points x_{1} such that, starting the Secant Method at x_{0} and any of these possible x_{1}’s, we reach a zero of f within one step.

46. i. Construct a function with a single zero in R such that for some (x_{0}, x_{1}) the Secant Method yields a 4-periodic orbit (namely (x_{4}, x_{5}) = (x_{0}, x_{1})).

ii. Construct a function, continuous on R and infinitely differentiable on R − {0}, having a single zero at the point 0, and possessing the following property: For any x_{0} ≠ 0 there exists a point x_{1} such that, starting the Secant Method at x_{0} and x_{1}, we obtain a 4-periodic orbit.

47. i. Find the general solution of the recursion d_{n+1} = 2d_{n} + C, where C is a constant.

ii. What does the first part hint regarding the error in Newton’s Method? iii. The same as the preceding two parts for the recursion d_{n+2} = d_{n} +

d_{n+1} + C and the Secant Method, respectively.

48. i. Suppose a new method for solving non-linear equations has been proposed. Given any three points, the process replaces the first of these by a new one, determined by the values of the given function at the initial points. The error in this method satisfies e_{n+3} ≈ Ke_{n+2}e_{n+1}e_{n}^{2} for some constant K. Compare the new method with the others we studied.

ii. Same for a method based on the last four points at each stage, and satisfying e_{n+4} ≈ Ke_{n+3}^{2}e_{n+2}e_{n+1}^{4}e_{n}^{6}.

Interpolation

49. Write down the Taylor polynomial of order n (around 0) of the function \frac{1}{x^{3}+1}. At what interval can it be expected to form a good approximation to the function?

50. i. Prove that there does not exist a polynomial P approximating the function e^{x} to within a fixed error on the whole line (namely, such that |e^{x} − P (x)| < C for every x ∈ R, where C is a constant).

ii. Prove that, if ε > 0 is sufficiently small, then there does not exist a polynomial P approximating the function sin x to within an error of ε on the whole line.

51. Let f be a continuous function on the interval [−a, a]. Use Weierstrass Theorem to prove the following:

i. If f is an even function (i.e., f(−x) = f(x) for every x), then for every ε > 0 there exists a polynomial P of the form P (x) = ∑_{i=0}^{n} c_{i}x^{2i}, satisfying |f(x) − P (x)| < ε for every x ∈ [−a, a].

ii. If f is odd, then there exist approximating polynomials of the form P (x) = ∑_{i=0}^{n} c_{i}x^{2i+1}.

52. Use Weierstrass Theorem to prove the following:

i. If f is a continuous function on [1, ∞) and f(x) ⟶_{x→∞} a, then for every ε > 0 there exists a polynomial P such that |f(x) − P (1/x)| < ε for every x ∈ [1, ∞).

ii. If f is a continuous function on [0, ∞) and f(x) ⟶_{x→∞} a, then for every ε > 0 there exists a polynomial P such that |f(x) − P (e^{−x})| < ε for every x ∈ [0, ∞).

53. Let f be a continuous function defined on [a, b] and (x_{k})_{k=1}^{∞} a sequence of distinct points in the interval. For each n, let P_{n} denote the interpolation polynomial of degree at most n−1 coinciding with f at the points x_{1}, x_{2}, …, x_{n}.

i. Prove that, if f is a polynomial, then P_{n} = f for sufficiently large n.

ii. Show by an example that, in general, the sequence (P_{n})_{n=1}^{∞} does not even need to converge to f as n ⟶ ∞.

54. Let f be a function defined on some interval [a, b] and x_{0}, x_{1}, …, x_{n} be distinct points of the interval, such that f is differentiable (at least) k_{i} times at x_{i} for 0 ≤ i ≤ n (where k_{i} ≥ 0 for each i). Prove that there exists a polynomial P having the following property:

P^{(k)}(x_{i}) = f^{(k)}(x_{i}), 1 ≤ i ≤ n, 0 ≤ k ≤ k_{i} .

What upper bound can you ensure for the degree of P?

55. Suppose f is defined on an interval [a, b] symmetric around 0 (i.e., a = −b), and the points x_{0}, x_{1}, …, x_{n} are symmetric around 0 (i.e., x_{n−i} = −x_{i} for each i). Show that, if f is an even (or odd) function, then so is the interpolation polynomial.

56. For each of the following functions determine for which values of n the Taylor polynomial of order n is the interpolation polynomial of the function for appropriately selected points:

i. f(x) = e^{x}, x ∈ [0, 2]. ii. f(x) = \frac{1}{1+x}, x ∈ [−\frac{1}{2}, \frac{1}{2}]. iii. f(x) = ln(1 + x), x ∈ [−\frac{1}{2}, \frac{1}{2}].

57. i. Find the polynomial P_{3}(x) interpolating the function (\frac{x}{π} + 1) cos x at the points 0, \frac{π}{4} , \frac{π}{3} and \frac{π}{2} .

ii. Compare the actual error obtained by predicting the value of the function at the point \frac{π}{6} using the interpolation polynomial with the bound on the error guaranteed by the theory.

58. For each of the following functions, check whether any sequence (P_{n})_{n=0}^{∞} of interpolation polynomials, interpolating the function at x_{0}^{(n)}, x_{1}^{(n)}, …, x_{n}^{(n)}, converges uniformly to the function:

i. f(x) = xe^{x}, x ∈ [a, b]. ii. f(x) = x^{2}e^{x}, x ∈ [a, b]. iii. f(x) = \frac{1}{1+x}, x ∈ [0, \frac{1}{2}]. iv. f(x) = \frac{1}{1−x}, x ∈ [0, 1).

59. For an arbitrary fixed sequence (x_{n})_{n=0}^{∞}, let P_{n} interpolate the function sin x at the points x_{0}, x_{1}, …, x_{n}.

i. Show that the sequence (x_{n})_{n=0}^{∞} can be constructed so that the sequence P_{n} converges to the function pointwise on the whole real line (that is, P_{n}(x) ⟶_{n→∞} sin x for every x ∈ R).

ii. Show that (x_{n})_{n=0}^{∞} can be selected so that P_{n} does not converge to the function pointwise on the whole line (i.e., there exists at least one point at which the convergence does not hold).

60. Same as the preceding question for the function e^{x}.

61. Same as the preceding question for the function e^{x^{2}}.

62. i. Use MATLAB to plot the first several interpolation polynomials, based on interpolating Runge’s function

f(x) = \frac{1}{1 + x^{2}} , −5 ≤ x ≤ 5,

at equally spaced points, and verify that they form poor approximations to the function.

ii. Do the same, interpolationg at Chebyshev’s points

x_{j}^{(n)} = 5 cos \frac{jπ}{n} , j = 0, 1, …, n .

Verify that this time the polynomials form good approximations to the function.

63. i. Prove the theorem stated in class regarding the error in the interpolation polynomial for equally spaced points in the cases n = 2 and n = 3.

ii. State and prove the analogue for n = 4.

64. Let P_{n} interpolate f at the points x_{0}, x_{1}, …, x_{n}. Find a simple expression for the polynomial interpolating the function xf(x) at the same points.

65. Let P_{n} interpolate f at the points x_{0}, x_{1}, …, x_{n}. Find a simple expression for the polynomial interpolating f at the points x_{1}, x_{2}, …, x_{n}.

66. Recalculate the polynomial P_{3} of Problem 57 using Newton’s divided differences.

67. Write a MATLAB function implementing Horner’s rule for evaluating the value of a polynomial at a given point.

68. i. Write a MATLAB function which, given a function f and interpolation points x_{0}, x_{1}, …, x_{n}, finds the Newton form of the interpolation polynomial. Minimize the time and space requirements.

ii. Write a MATLAB function which, given a polynomial P in Newton form, finds its usual representation. Minimize time and space.

iii. Write a MATLAB function which, given a polynomial P in Newton form and a point t, calculates the value of P at t. Minimize time and space.

69. Let n be an arbitrary fixed non-negative integer, and t_{0}, t_{1}, …, t_{n} be all zeros of Chebyshev’s polynomial T_{n} of degree n.

i. Prove the following properties of Chebyshev’s polynomials (T_{i})_{i=0}^{n}: a. ∑_{k=0}^{n} T_{i}(t_{k})T_{j}(t_{k}) = 0, i ≠ j .

b. ∑_{k=0}^{n} T_{i}^{2}(t_{k}) = \{ \begin{matrix}n + 1, & i = 0, \\ \frac{n+1}{2}, & i = 1, 2, …, n . \end{matrix}

ii. Interpret the preceding part as providing a certain orthonormal basis of R^{n+1}.

70. Use the techniques of recurrence sequences to obtain an “explicit” expression for T_{n}(x).

71. Find all zeros of T′_{n}, and use them to obtain an explicit expression for T′_{n}(x) as a product of linear factors.

Curve Fitting

72. Prove that the error measures presented in class for testing the compatibility of a function to given data always satisfy E_{1}(f) ≤ E_{2}(f) ≤ E_{∞}(f), with equalities if and only if f(x_{1}) − y_{1} = f(x_{2}) − y_{2} = … = f(x_{n}) − y_{n}.

73. Find the least-squares line for the following data points:

i. (−1, 2), (0, 3), (1, 5), (2, 4) . ii. {(k, k^{2}) : k = 0, 1, …, n} .

74. Prove that:

i. If the data points are symmetric with respect to the origin, then the least-squares line passes through the origin.

ii. If the data points are symmetric with respect to the y-axis, then the least-squares line is horizontal.

75. Prove that x̄^{2} ≥ x̄^{2}, with equality if and only if x_{1} = x_{2} = … = x_{n}.

76. Let y = ax + b be the least-squares line corresponding to the points (x_{1}, y_{1}), (x_{2}, y_{2}), …, (x_{n}, y_{n}). Given constants c_{1}, c_{2}, d_{1}, d_{2}, find the least-squares line corresponding to the points (c_{1}x_{k} + d_{1}, c_{2}y_{k} + d_{2}), k = 1, 2, …, n. Explain your result intuitively.

77. Given data points (x_{1}, y_{1}), (x_{2}, y_{2}), …, (x_{n}, y_{n}), find the equations defining the least-squares fit for each of the following families of curves:

i. Lines through the origin y = ax. ii. Parabolas of the form y = ax^{2}. iii. Parabolas of the form y = ax^{2} + b.

iv. The family of curves y = a sin x + b cos x.

78. Consider the example given in class of the least-squares line corresponding to the points (0, 0), (1, 1), (2, 4), (3, 9). Suppose the y-coordinates of these points are perturbed in such a way that the sum total of all changes does not exceed 1. That is, if the new y-coordinates are y_{0} through y_{3}, then |y_{0} − 0| + |y_{1} − 1| + |y_{2} − 4| + |y_{3} − 9| ≤ 1.

i. Find the range of possible slopes a of the corresponding least-squares lines y = ax + b.

ii. Same for the range of possible values of b. iii. Find the set of all possible values of the pair (a, b).

79. For each of the following families of curves use data linearization to find a general formula for the fit corresponding to given data points:

i. y = \frac{a}{x+b}.

ii. y = ax^{b}. iii. y = \frac{1}{(ax+b)^{2}} .

iv. y = \frac{1}{ax^{2}+b}.

v. y = \frac{1}{1+be^{ax}}.

80. For each of the following families of curves, write down the conditions on the parameters necessary for the curve to be the least-squares fit. Obtain explicit expressions wherever possible, or at least try to simplify the equations you got:

i. y = \frac{a}{x+b}. ii. y = x^{a} + b. iii. y = a sin x + b. iv. y = a + be^{x} + ce^{−x}. v. y = (x + a)^{2}. vi. y = a sin^{2} x + b cos^{2} x.

81. The least-squares plane z = ax + by + c corresponding to the n data points (x_{1}, y_{1}, z_{1}), …, (x_{n}, y_{n}, z_{n}) is the plane obtained by minimizing

E_{2}(a, b, c) = (\frac{1}{n} ∑_{k=1}^{n}(ax_{k} + by_{k} + c − z_{k})^{2})^{1/2} .

i. Find the equations determining the parametrs of the least-squares plane in the general case.

ii. Find the least-squares plane corresponding to the data points (1, 4, 5), (2, 7, 9), (3, 9, 12), and (4, 12, 10).

iii. Find the least-squares plane for the data points {(k, k^{2}, k^{3}) : k = 0, 1, …, n} .

82. Consider the system of equations for the parameters m_{i}, i = 0, 1, …, n, required for specifying the cubic spline, in each of the following cases, and show that in each of them the matrix of coefficients is diagonally-dominant.

i. The first derivative at the endpoints is known: m_{0} = \frac{3}{h_{0}} (d_{0} − S′(x_{0})) − \frac{m_{1}}{2} , m_{n} = \frac{3}{h_{n−1}} (S′(x_{n}) − d_{n−1}) − \frac{m_{n−1}}{2} .

ii. The second derivative is extrapolated to the endpoints: m_{0} = m_{1} − \frac{h_{0}}{h_{1}} (m_{2} − m_{1}), m_{n} = m_{n−1} + \frac{h_{n−1}}{h_{n−2}} (m_{n−1} − m_{n−2} .

iii. The second derivative is assumed to be nearly constant near the endpoints:

m_{0} = m_{1}, m_{n} = m_{n−1} .

iv. The second derivative is assumed to vanish at the endpoints: m_{0} = m_{n} = 0 .

83. Find the cubic splines passing through the points (1,1), (2,3), (3,9) and (4,7), obtained for each of the boundary conditions specified in the preceding question.

84. Suppose instead of piecewise-cubic twice continuously differentiable splines we look for piecewise-quartic three times continuously differentiable splines. How many degrees of freedom should be expected when selecting such splines?

85. Prove that any diagonally-dominant matrix A ∈ M_{n}(C) is non-singular, following the indicated steps:

i. Multiplying A by an appropriately chosen matrix, show that, without loss of generality, we may assume all elements on the main diagonal of A to be equal to 1. Thus A = I + T, where I is the identity matrix and T = (t_{ij})_{i,j=1}^{n} satisfies ∑_{j=1}^{n} |t_{ij}| < 1 for 1 ≤ i ≤ n. (In fact, the diagonal entries of T are 0, but this is of no consequence here.)

ii. For D = (d_{ij})_{i,j=1}^{n} in M_{n}(C), denote ∥D∥ = max_{1≤i≤n} ∑_{j=1}^{n} |d_{ij}| . Prove the following properties of ∥ · ∥:

a. ∥D∥ ≥ 0 for every D ∈ M_{n}(C), with equality if and only if

D = 0.

b. ∥αD∥ = |α| · ∥D∥ for every α ∈ C and D ∈ M_{n}(C).
c. ∥D_{1} + D_{2}∥ ≤ ∥D_{1}∥ + ∥D_{2}∥ for every D_{1}, D_{2} ∈ M_{n}(C).
d. ∥D_{1}D_{2}∥ ≤ ∥D_{1}∥ · ∥D_{2}∥ for every D_{1}, D_{2} ∈ M_{n}(C).

iii. A sequence (D_{m})_{m=1}^{∞} of matrices in M_{n}(C) converges to D ∈ M_{n}(C), and we denote D_{m} ⟶_{m→∞} D, if ∥D_{m} − D∥ ⟶_{m→∞} 0 as a sequence of real numbers. Prove that D_{m} ⟶_{m→∞} D if and only if we have entrywise convergence.

iv. Prove that, if D_{m} ⟶_{m→∞} D and E_{m} ⟶_{m→∞} E, then D_{m}+E_{m} ⟶_{m→∞} D+E

and D_{m}E_{m} ⟶_{m→∞} DE.

v. Prove that, if D ∈ M_{n}(C) with ∥D∥ < 1, then the sequence of matrices (I + D + D^{2} + … + D^{m−1})_{m=1}^{∞} converges.

vi. Conclude the proof by showing that A^{−1} = I − T + T^{2} − T^{3} + ….

86. Write a MATLAB function which, given points (x_{0}, y_{0}), (x_{1}, y_{1}), … (x_{n}, y_{n}) and required values of the first derivative at the endpoints, finds the corresponding cubic spline.

87. Test the theorem regarding the minimality property of cubic splines experimentally as follows. Fix n (say n = 10). Find the integral of the square of the second derivative of both the interpolation polynomial and the cubic spline having the same first derivative as the interpolation polynomial at the endpoints of the interval for many (say 1,000) sets of (n+1) points {(x_{k}, y_{k}) : k = 0, 1, … n} in the plane, constructed according to the following mechanism:

a. x_{0} = y_{0} = 0.

b. Suppose (x_{k}, y_{k}) has been selected. Choose randomly and independently two numbers, H_{1} in the interval [0, 1] and H_{2} in the interval [−1, 1]. Put x_{k+1} = x_{k} + H_{1} and y_{k+1} = y_{k} + H_{1}H_{2}.

Numerical Differentiation

88. Suppose f is defined only to the right of the point x, so that the central difference approximation for f′(x), namely

f′(x) = \frac{f(x + h) − f(x − h)}{2h} + O(h^{2}) ,

is inapplicable. Can you suggest an alternative O(h^{2}) approximation for the required (one-sided) derivative?

89. Prove the following approximation for f^{(4)}(x), mentioned in class:

f^{(4)}(x) = \frac{−f_{3} + 12f_{2} − 39f_{1} + 56f_{0} − 39f_{−1} + 12f_{−2} − f_{−3}}{6h^{4}} + O(h^{4}) ,

where f_{i} = f(x + ih) for each i.

90. Prove that for any fixed non-zero distinct numbers λ_{1}, λ_{2}, …, λ_{k}, one can find an O(h^{k}) approximation for f′(x) in terms of f(x), f(x + λ_{1}h), f(x + λ_{2}h), …, f(x + λ_{k}h).

Numerical Integration

91. Show that, viewing the constant function f(a) as an interpolation polynomial of degree 0 for the function f(x), and using the formula for the interpolation error, we can arrive at the bound on the error in the rectangle method obtained in class.

92. The integral ∫_{0}^{1} \frac{e^{x}}{x+1} is to be evaluated to within an error of 0.001 by the rectangle method. Into how many subintervals n do you need to divide [0, 1] to obtain the required precision? (Try to minimize n. Ignore computer errors.)

93. How can you approximate the integral ∫_{0}^{1}(1 − \sqrt{2x − x^{2}})dx to within an error of 0.001 by the rectangle method?

94. Consider the function log on the interval [\frac{1}{2}, n + \frac{1}{2}].

i. Find the exact value of the integral ∫_{1/2}^{n+1/2} logxdx.

ii. Estimate the integral employing the midpoint rule, with division to n subintervals. Using the estimate for the error, obtain a weak form of Stirling’s formula (providing an asymptotic expression for n!).

iii. What can you get in the previous part if you use the bound on the error for each subinterval separately, and not just for the whole interval?

iv. Repeat the above for the integral ∫_{1/2n}^{1+1/2n} logxdx.
95. It is well known that, denoting H_{n} = 1 + \frac{1}{2} + \frac{1}{3} + … + \frac{1}{n}, we have

lim_{n→∞}(H_{n} − logn) = γ,

where γ is Euler’s constant. How much of this do you obtain by approximating the integral ∫_{1/2}^{n+1/2} \frac{dx}{x} using the midpoint rule, with division to n intervals, and estimating the error?

96. Same as Problem 91 for the interpolation polynomial of degree 1 based on the values of the function at both endpoints of the interval and the trapezoidal rule.

97. Verify directly that, using the trapezoidal rule with division to n subintervals to estimate the integral ∫_{a}^{b} x^{2}dx, the resulting error is within the predicted bounds.

98. Verify directly that, using Simpson’s rule with division to n subintervals to estimate the integral ∫_{0}^{1} e^{x}dx, the resulting error is within the predicted bounds.

99. i. Prove that there exists no approximation formula for ∫_{a}^{b} f(x)dx, based on the values of f at finitely many points only, which is completely precise for all polynomials.

ii. Prove that, moreover, for any such approximation formula, based on the values of f at n points, there exists a polynomial of degree at most 2n for which the formula gives an imprecise answer.

iii. Explain why the bound 2n in the preceding part is best possible.

100. Show that Boole’s rule gives the exact value for integrals of polynomials of degree not exceeding 5.

101. Let (P_{n}(x))_{n=0}^{∞} be a sequence of polynomials. Show that if deg P_{n} = n for each n then {P_{n}(x) : n = 0, 1, 2, …} is a basis of the space of all polynomials.

102. Verify the following properties of Legendre polynomials:

i. ϕ′_{n}(x) = xϕ′_{n−1}(x) + nϕ_{n−1}(x).

ii. ϕ_{n}(x) = xϕ_{n−1}(x) + \frac{x^{2}−1}{n}ϕ′_{n−1}(x)).

iii. (n + 1)ϕ_{n+1}(x) = (2n + 1)xϕ_{n}(x) − nϕ_{n−1}(x).

iv. ∫_{−1}^{1} ϕ_{n}(x)^{2}dx = \frac{2n−1}{2n+1} ∫_{−1}^{1} ϕ_{n−1}(x)^{2}dx.

v. ∫_{−1}^{1} ϕ_{n}(x)^{2}dx = \frac{2}{2n+1}.

vi. The set of zeros of ϕ_{n} is symmetric around 0.

103. Carry out the Gram-Schmidt process for the system {1, x, x^{2}, …} on the interval [−1, 1] to find the polynomials ϕ_{n}(x), 0 ≤ n ≤ 4.

103. i. Use Gaussian integration with two and with three division points to estimate ∫_{−1}^{1} \frac{dx}{1+x^{2}} . Compare with the correct value and verify, for two division points, that the error is within the allowed limits.

ii. Same for the integral ∫_{−1}^{1} \frac{dx}{x+2}.

104. Use Gaussian integration with two and with three division points to estimate ∫_{−1}^{1} x^{n}dx. What are the absolute and the relative errors? Explain!

105. Use Gaussian integration with two division points to estimate the integral ∫_{−1}^{1} cos πnx dx. Use MATLAB to check how good the estimate is for various values of n. Do you think that the points ±\frac{1}{\sqrt{3}} used for the estimate are more representative for these integrals (for large n) than random points in [−1, 1]? Explain!

106. Let F denote the vector space of all Riemann-integrable functions over a finite interval [a, b], and let w ∈ F be strictly positive throughout the interval.

Show that the function ⟨·, ·⟩ : F × F ⟶ R defined by

⟨f, g⟩ = ∫_{a}^{b} f(x)g(x)w(x)dx, f, g ∈ F

forms an inner product on F.

107. Verify that each of the following defines an inner product on the space of polynomials:

a. ⟨P, Q⟩ = ∫_{0}^{∞} P (x)Q(x)e^{−x}dx.

b. ⟨P, Q⟩ = ∫_{−∞}^{∞} P (x)Q(x)e^{−x^{2}}dx.
c. ⟨P, Q⟩ = ∫_{−1}^{1} \frac{P (x)Q(x)}{\sqrt{1−x^{2}}}dx.

108. i. Using the equality ∫_{−∞}^{∞} e^{−x^{2}}dx = \sqrt{π} show that one can calculate ∫_{−∞}^{∞} e^{−x^{2}}P (x)dx precisely for any polynomial P.

ii. Calculate ∫_{0}^{∞} \sqrt{x}e^{−x}dx.

109. Suppose we are interested in integrals of the form ∫_{0}^{1} f(x) ln \frac{1}{x}dx.

i. Write down the equations obtained by trying to find an approximation of the form

I = w_{1}f(x_{1}) + w_{2}f(x_{2}) + … w_{k}f(x_{k})

for this integral, which will be precise if f is a polynomial of degree not exceeding 2k − 1.

ii. Define an appropriate inner product on the space of continuous realvalued functions on [0, 1], which will be of use for finding the division points x_{1}, x_{2}, …, x_{k}.

iii. Find the division points for k = 2.

Review Problems

1. i. A system employs the IEEE standard with double precision, round-to-nearest mode and correctly rounded arithmetic. Find two floating point numbers a, b ∈ [0, 1] for which the rounding error when computing a + b is as large as possible.

ii. System as before, but with round-towards-∞. Find two floating point numbers a, b ∈ [0, 1] for which the rounding error when computing a·b is as large aspossible.

2. i. Find the rate of convergence of Newton’s Method for each root nπ of the function f(x) = x sin x.

ii. For each such root, find explicitly an interval [nπ − ε_{1}, nπ + ε_{2}] such that, starting Newton’s Method at any point of the interval, we obtain a sequence converging to nπ.

3. Let (x_{k})_{k=1}^{∞} be a sequence of distinct points on the real line. For each n, let P_{n} denote the interpolation polynomial of degree at most n − 1 coinciding with the function f(x) = x + sin x at the points x_{1}, x_{2}, …, x_{n}.

i. Prove that we do not necessarily have P_{n}(x) ⟶_{n→∞} f(x) for every x ∈ R.

ii. Prove that, if the sequence x_{k} is bounded, then P_{n}(x) ⟶_{n→∞} f(x) for every x ∈ R, and, moreover, the convergence is uniform on any finite interval [a, b].

iii. Does there exist any sequence of polynomials converging to f uniformly on R?

4. The interpolation polynomial P of degree at most n satisfying P (x_{i}) = f_{i} for 0 ≤ i ≤ n has been calculated. It turned out that, in fact, the value at some x_{i} (for a single i) is not f_{i}, but rather another number v. Let P_{1} denote the updated interpolation polynomial. Given a point a, explain how you can calculate P_{1}(a) as fast as possible. How much time is required for this calculation (up to a constant factor)?

5. Let f be a twice differentiable function defined on the interval [0, 2]. It is known that |f′′(x)| ≤ M_{1} for x ∈ [0, 1] and |f′′(x)| ≤ M_{2} for x ∈ [1, 2], where M_{1} and M_{2} are given positive constants. You are required to estimate ∫_{0}^{2} f(x)dx using the trapezoidal rule, dividing each of the two intervals [0, 1] and [1, 2] into equal subintervals. If you are allowed to calculate f at n + 1 points altogether, into how many subintervals will you subdivide [0, 1] and into how many [1, 2]?