Coursework [100 marks]

Question 1.  [20 marks]  Self-consistent iteration.

Solve the transcendental equation

x = e −x

using fixed-point iteration. That is, using the initial guess x0  = 1, obtain a sequence of real numbers

x1  = e −x0

x2  = e −x1 

x50  = e −x49 

which tends to the value x∞  = 0.567143 ..., that is, a root of the function f(x) = x − e −x . Formally, the iteration can be written as

xi+1  = e −xi  for i = 0, 1, 2, . . . with x0  = 1.

The limit x∞ satisfies f(x∞ ) = 0.

(a)  Write a for loop that performs the above iteration, starting from the initial condition    x0  = 1. Use an if and a break statement to stop the loop when the absolute value of the difference |xi+1 − xi | between two consecutive iterations is less than a prescribed    tolerance ! = 10 − 15 . Use a long double to store the values of x, and the change in  x, between iterations. Use setprecision(18) to print out the final value of x to 18

digits of accuracy.                                                                                                                     [10] (b)  In how many iterations did your loop converge?                                                                      [5]

(c)  What is the final error in the above transcendental equation? [Hint: use the final value of       x to compute and display the difference x − e −x .] Is the error what you expected (i.e. is        it smaller than !)?                                                                                                                        [5]

Question 2.  [20 marks] Inner products, sums and norms.

(a)  Write a function named dot that takes as input two vectors A = {A1 , ...,An } ∈ Rn and B = {B1 , ...,Bn } ∈ Rn (of type valarray<long double>) and returns their inner

product

n

A · B =  Ai Bi

as a long double number.

(b)  Use this function to compute the sum       for n = 106 . To do so, you may define a

Euclidean n-vector A={1,1/2,1/3, . . .,1/1000000} as a valarray<long

double> and compute the inner product A · A. Print the difference A · A − π  /62

between your result and the exact value of the infinite sum on the screen.       Remark: Define π = 3.1415926535897932385 to long double accuracy.

(c)  Repeat Question 2(a) using compensated (Kahan) summation and name your function cdot. Use the old dot function and the new cdot function to compute the sum          c2 for n = 106 and c = 0.1.

Hint: Define a constant Euclidean n-vector C={c,c,...,c} as a valarray<long double> of size n and compute the inner product C · C.

Print the difference between your result and the exact value nc2 .                                           [5]

(d)  Write code for a function object that has a member variable m of type int, a suitable constructor, and a member function of the form

double operator()(const valarray<long double>  A) const { which returns the weighted norm

lm (A) =  m%#$$ |Ai |m  = &  |Ai |m '1/m

(2)

Use this function object to calculate the norm l2 (A) for the n-vector A in Question 2b.           Does the quantity l2 (A)2 equal the inner product A · A that you obtained above? [Note:         half marks awarded if you use a regular function instead of a function object.]                    [5]

Question 3.  [20 marks] Finite differences.

(a)  Write a C++ program that uses finite difference methods to numerically evaluate the

second derivative f (x) of a function f(x) whose values on a fixed grid of points are     specified f(xi ) = fi , i = 0, 1, ...,N . Your code should use valarray<long      double> containers to store the values of xi , fi and f . Assume the grid-points are     located at xi  = (2i − N)/N on the interval [−1, 1] and use 2nd order finite differencing to compute an approximation f for f (xi ):

f   = f   =

f   =

fi+2 − 2fi+1 + fi

∆x2

fi+1  − 2fi + fi − 1

∆x2

fi  − 2fi − 1 + fi −2

∆x2

if   i = 0

if   i = 1, 2, ...,N − 1

if   i = N

with ∆x = 2/N. Demonstrate that your program works by evaluating the second         derivatives of a known function, f(x) = e − 16x2 , with N +1 = 64 points. Compute the difference between your numerical derivatives and the known analytical ones:

ei  = fumerical (xi ) − f (xi )

double> on the screen and tabulate (or plot) them in your report to 10 digits of

accuracy. Comment on whether the error is what you expected.                                            [10]

(b)  For the same choice of f(x), demonstrate 1st-order convergence, by showing that, as N

increases, the mean error  e$ decreases proportionally to ∆x = 2/N . You may do so by tabulating the quantity N e$ for different values of N (e.g. N +1 = 16, 32, 64, 128)     and checking if this quantity is approximately constant. (Alternatively, you may plot      log e$ vs. log N and check if the dependence is linear and if the slope is −1.)

Here, the mean error  e$ is defined by

N

e$ = ! |ei | = l1 ( ).

Hint: you may use your code from Question 2d to compute the mean error.                       [10]

Question 4.  [20 marks] Numerical integration.

We wish to compute the definite integral

I =  b       (b − x)xdx

numerically, with endpoints a = 0 and b = 4, and compare to the exact result, Iexact  = 2π . In Questions 4a, 4b and 4c below, use instances of a valarray<long double> to store the values of the gridpoints xi , function values fi  = f(xi ) and numerical integration weights wi .

(a)  Use the composite trapezium rule

ab f(x)dx   wi fi  = w · f ,    w =  ∗ {1, 2, 2, ..., 2, 1}

to compute the integral I, using N +1 = 64 equidistant points in x ∈ [a,b], that is,

xi  = a + i∆x, for i = 0, 1, ...,N, where ∆x = ba  is the grid spacing. Output your               numerical result Itrapezium  and the difference Itrapezium  − Iexact . [Hint: Use the function              dot or cdot, created in Question 2a or 2b, to easily evaluate the inner product w · f].     [5]

(b)  Use the extended Simpson rule

ab f(x)dx   wi fi ,    w =  ∗ {17, 59, 43, 49, 48, 48, ..., 48, 49, 43, 59, 17}

your numerical result ISimpson  and the difference ISimpson  − Iexact .                                              [5]

(c)  Use the Clenshaw-Curtis quadrature rule

ab f(x)dx   wi fi ,

wi  = ∗* 2 1 −"1)/2

 ,

i = 0 or i = N 1 ≤ i ≤ N − 1

,

on a grid of N +1 = 64 points xi  = [(a + b)+(a − b)cosθi]/2, where θi  = iπ/N, i = 0, 1, ...,N to compute the integral I.

Hint: First compute the values of θi , xi , fi  = f(xi ) and wi and store them in                  valarray<long double> containers, as shown in the lab. Then, you may use the function from Question 2a or 2c to compute the inner product w · f .

Output to the screen (and list in your report) your numerical result IClenshawCurtis  and the difference IClenshawCurtis  − Iexact .

(d)  Compute the integral I using a Mean Value Monte Carlo method with N = 10000     samples. Output to the screen (and list in your report) your numerical result IMonteCarlo and the difference IMonteCarlo  − Iexact .

[5]

[5]