Coursework [100 marks]

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

Using a pocket calculator, one may notice that by applying the cosine key repeatedly to the value (in radians) one obtains a sequence of real numbers

x1  = cosx0  = 1.

x2  = cosx1  = 0.54030230586814

x3  = cosx2  = 0.54030230586814

x4  = cosx3  = 0.857553215846393 

x21  = cosx20  = 0.739184399771494 

which tends to the value x∞  = 0.739085 ..., which is the point where the graphs of the function x and cos x intersect. The iteration can be written as

xn+1  = cosxn for n = 0, 1, 2, . . . with x0  = 1.

The limit x∞ satisfies the transcendental equation

cos x = x.

Write a for loop that performs the iteration xn+1  = cosxn starting from the initial condition x0  = 0 and stops when the absolute value of the difference |xn+1 − xn | between two                consecutive iterations is less than a prescribed tolerance e = 10 −12 . Use an If/Break          statement to exit the loop when the condition is met. Print out the final value xn+1 to 16 digits of accuracy. In how many iterations did your loop converge? What is the final error in the       above transcedental equation? (Hint: use the final value to compute and print out the

difference x − cosx.)                                                                                                                         [20]

Question 2.  [20 marks]  Inner products.

(a)  Write a function that takes as input two Euclidean vectors ⃗u = {u1 ,u2 , . . . ,uN } ∈ RN and ⃗v = {v1 ,v2 , . . . ,vN } ∈ RN  (of type valarray<long  double>) and returns  their inner product (also known as Hadamard product)

N

⃗u · ⃗v =工 uivi                                                                            (1)

i=1

as a long  double number. Your function may use  (u *v) .sum() to compute the dot product of the valarrays u and v. Create a constant valarray<long           double> equal to u={0 . 1,0 . 1, . . . ,0 . 1} with N = 106 elements. Demonstrate

Display the difference ⃗u · ⃗u − 104 on the screen.                                                                     [5]

(b)  Repeat Question 2a using Kahan compensated summation to compute the sum.                  [5]

(c)  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<double>   u)  const  {

}    ...

which returns the weighted norm

「

' N

ℓm (⃗v) =  m'/ |vi |m

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

Question 3.  [20 marks] Finite differences.

(a)  Write a C++ program that uses finite difference methods to numerically evaluate the       first derivative of a function f(x) whose values on a fixed grid of points are specified      f(xi ),  i = 0, 1, . . . ,N . Your code should use three instances of a valarray<long      double> to store the values of xi , f(xi ) and f\ (xi ). Assume the grid-points are located at xi  = a + i∆x with ∆x = (b − a)/N . on the interval x ∈ [a,b] and use 2nd order       finite differencing to compute an approximation for f\ (xi ):

−3f(x0 ) + 4f(x1 ) − f(x2 )

f\ (xi )   =    + O(∆x2 )   for    i = 1, 2, . . . ,N − 1

f(xN−2) − 4f(xN−1) + 3f(xN )

2∆x

Demonstrate that your program works by evaluating the derivatives of a known function, f(x) = sin3x, with N + 1 = 32 points on the interval x ∈ [a,b] = [−1, 1]. Compute the difference between your numerical derivatives and the known analytical ones:

ei  = fumerical (xi ) − fnalytical (xi )

at each grid-point. Output the values ei of this valarray<long  double> on the             screen and tabulate (or plot) them in your report.                                                                  [10]

(b)  For the same choice of f(x), demonstrate 2nd-order convergence, by showing that, as N increases, the mean error ⟨e⟩ decreases proportionally to ∆x2  ∝  N −2  . You may do so  by tabulating the quantity N2 ⟨e⟩ for different values of N (e.g. N + 1 = 16, 32,

64,  128) and checking if this quantity is roughly constant. Alternatively, you may plot log⟨e⟩ vs. log N and check if the dependence is linear and if the slope is -2.

Here, the mean error ⟨e⟩ is defined by

N

⟨e⟩ =  |ei | = ℓ 1 (⃗e).

You may use your code from Question 3c to compute the ℓ1 norm.                                     [10]

Question 4.  [20 marks]  Numerical integration.

We wish to compute the definite integral

I = \a b sin ( ) dx

numerically for a = 0, b = 10 and compare to the exact result, Iexact  = 2.74324739415100920.

(a)  Use the composite trapezium rule

\a b f(x)dx ≃ wi fi ,    wi  =     ∆x =  ,

to compute the integral I, using N + 1 = 128 equidistant points in x ∈ [a,b]. Use three instances of a valarray<long  double> to store the values of the gridpoints xi ,    function values fi  = f(xi ) and weights wi . [Hint: you may use the function from          Question 2a to compute the dot product of the valarrays wi and fi .] Output to the    screen (and list in your report) your numerical result Itrapezium  and the difference

Itrapezium  − Iexact .                                                                                                                          [5]

(b)  Use the composite Hermite rule

\a b f(x)dx ≃  wi fi +  [f\ (a) − f\ (b)]

with the derivatives f\ (x) at x = a and x = b evaluated analytically (and the weights wi identical to those given above for the trapezium rule), to compute the integral I, using   N + 1 = 128 equidistant points in x ∈ [0, 1]. Output to the screen (and list in your

report) your numerical result IHermite  and the difference IHermite  − Iexact .                                [5]

(c)  Use the Clenshaw-Curtis quadrature rule

\a b f(x)dx ≃  wi fi ,

wi  = ∗ { 2 ,(1 −z1(−)1)/2

i = 0 or i = N

)   1 ≤ i ≤ N − 1 ,

on a grid of N + 1 = 128 points xi  = [(a + b) − (b − a)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 as valarrays. Then, you may use the function from Question 2a to compute the dot product of the valarrays wi and fi .] Output to the     screen (and list in your report) your numerical result IClenshawCurtis  and the difference

(d)  Compute the integral I using a Mean Value Monte Carlo method with N = 1000,

N = 10000 and N = 100000 samples. Output to the screen (and list in your report)

your numerical results IMonteCarlo  and the difference IMonteCarlo  − Iexact for each N.               [5]

Question 5.  [20 marks]  Stellar equilibrium.

Consider the generalized Lane-Emden equation, in the form of the initial value problem

2x,                                                          (3)

This equation describes the structure of a self-gravitating planet or stellar object with mass density ρ(h) = h, where h is the specific enthalpy. Setting m(x) := −x2 h\ (x), the above   2nd-order differential equation can be reduced to a system of two 1st-order differential      equations for h(x) and m(x):

 h\ (x) =

〈 m\ (x) = h(x) x2

  h(0) = 1 

(m(0) = 0

The above equation appears singular at x = 0, but one can use a Taylor expansion about the origin to show that m(x) ≃ x3 + O(x5 ) for small x, so that h\ (0) = 0.

Declare a function

valarray<long  double>  F(const  long  double  t,  const

valarray<long  double>&  u)  {

}     ...

that takes the radius x and a valarray with elements ⃗u = {h,m} as arguments and          returns a valarray with components f⃗= {h\ ,m\ } as output, with h\ ,m\ given by Eq. (4)

above.                                                                                                                                           [5]

(b)  Solve the above 1st-order system numerically, with a 4th-order Runge-Kutta method

long  double  RK4(const  double  t,  const  double  dt,  const        valarray<long  double>   &,  long  double  f(const  long  double,

valarray<long  double>   &)  {

}     ...

using N + 1 = 151 equidistant points in x ∈ [0,π].

Declare a valarray of valarrays:

valarray<valarray<long  double>>  u

and use it to store ⃗u(xi ) = {hi ,mi } = {h(xi ),m(xi )} at each grid-point xi .

Output (only) the values {x0 ,x10 ,x20 , . . . ,xN } and {h(x0 ),h(x10 ),h(x20 ), . . . ,h(xN )} to         the screen and tabulate them in your report.                                                                             [5]

(c)  Compute the difference e(x) = hnumerical (x) − hexact (x) between your numerical solution hnumerical (x) and the exact solution

hexact (x) = sinc x = {

x  0

x = 0