Main Examination period 2019/20
MTH6150: Numerical Computing in C and C++ Assignment date: Thursday 23/4/2020
Submission deadline: Monday 8/6/2020 at 23:55
The coursework is due by Monday, 8th June 2020 at 23:55 BST. Please submit a report (in pdf format) containing answers to all questions, complete with written explanations and printed tables or figures. Tables should contain a label for each column. Figures should contain a title, axis labels and legend. Please provide a brief explanation of any original algorithm used in the solution of the questions (if not discussed in lectures). Code com- ments should be used to briefly explain what your code is doing. You need to show that your program works using suitable examples. Build and run your code with any of the free IDEs discussed in class (such as Visual Studio, Xcode, CLion, or the online compiler https://coliru.stacked-crooked.com/). The code produced to answer each question should be submitted aside with the report, in a single cpp file, or multiple cpp files compressed into a single zip file. It is useful to organise the code in different directo- ries, one for each question. The code for each question should also be copied into your report, so that it can be commented on when grading. This coursework should not take more than a total of 24 hours (assuming you have been engaging with lectures, labs and homework). Only material submitted through QMPlus will be accepted. Late submissions will be treated in accordance with current College regulations. Plagiarism is an assessment offence and carries severe penalties. See the accompanying guidelines for additional information.
Coursework [100 marks]
Question 1. [10 marks] Self-consistent iteration.
The geometry outside a spherical black hole can be expressed using a radial Schwarzschild coordinate r ∈ [1, ∞) or a radial tortoise coordinate x ∈ (−∞, ∞). The transformation between the two coordinate systems, x = x(r), is given by
x = r + ln |r − 1| . (1)
We seek the inverse transformation, r = r(x). For a given coordinate x, we wish to obtain the coordinate r satisfying Eq. (1). Write a C++ function rofx that takes x as argument, solves Eq. (1) for r via self-consistent iteration, and returns the value of r, up to a specified tolerance ε. The iteration can be performed as a for or while loop that computes the sequence
rn+1 =x−ln|rn −1| for n=0,1,2,… (2)
and stops when |rn+1 − rn| < ε. The function rofx should take x as input and may use
r0 = 0.8x as initial guess to perform the above iteration, in order to find r1, r2, r3, .... The
value of x should be fixed throughout the iteration, but the value of r should be updated until it converges. The function rofx should return the converged value of r. You may use a
tolerance parameter ε = 10−12. Run your code for a value of x > 2. In how many iterations
did it converge? What is the final value of r? Substitute r into Eq. (1): is it satisfied? [10]
1
Question 2. [10 marks] Inner products.
(a) Write a function that takes as input two vectors ⃗u = {u0, u1, …, uN } ∈ RN+1 and ⃗v = {v0,v1,…,vN} ∈ RN+1 (of typevector
XN
⃗u ∙ ⃗v =
as a double number. Demonstrate that your program works by computing the inner
u i v i
product of two real Euclidean 3-vectors, u={2,7,2} and v={3,1,4}.
( 3 )
(b) 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 vector
which returns the weighted norm
i=0
N
[5]
lm(⃗u,⃗v)= X|uivi| i=0
Question 3. [15 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 vector
xi = (2i − N )/N on the interval [−1, 1] and use 2nd order finite differencing to compute an approximation for f0(xi):
f0(x0) = −3f(x0) + 4f(x1) − f(x2) + O(Δx2) for i = 0 2Δx
f0(xi) = f(xi+1)−f(xi−1)+O(Δx2) for i=1,2,…,N−1 2Δx
f0(xN) = f(xN−2)−4f(xN−1)+3f(xN)+O(Δx2) for i=N 2Δx
with Δx = 2/N . Demonstrate that your program works by evaluating the derivatives of a known function, f (x) = e−x2 , with N + 1 = 16 points. Compute the difference between your numerical derivatives and the known analytical ones:
ei = fn0umerical(xi) − fa0nalytical(xi)
at each grid-point. Output the values ei of this vector
[8]
m
m/2
(4)
Use this function object to calculate the norms l1(⃗u,⃗v) and l2(⃗u,⃗v) for the 3-vectors in Question 2a. Does the quantity l2(⃗u, ⃗v)2 equal the inner product ⃗u ∙ ⃗v that you obtained above? [5]
2
(b) For the same choice of f (x), demonstrate 2nd-order convergence, by showing that, as N increases, the mean error ⟨e⟩ or the root mean square error ⟨e2⟩1/2 decrease proportionally to Δx2 ∝ N−2 . You may do so by tabulating (or log-log plotting) the quantity N2⟨e⟩ (or N2⟨e2⟩1/2) for different values of N (e.g. N + 1 = 16, 32, 64, 128). Is the quantity N2⟨e⟩ (or N2⟨e2⟩1/2) approximately constant?
Here, the mean error ⟨e⟩ is defined by
1 XN
⟨e⟩ = N + 1 i=0 |ei|. Alternatively, you may use the rms (root mean square) error
1 XN ⟨e2⟩1/2 = N + 1 i=0 |ei|2.
Hint: Given the vector ⃗e, you may use your C++ implementation of Eq. (4) from
Question 2 to compute ⟨e⟩ = l1(⃗e,⃗e)/(N + 1) or ⟨e2⟩1/2 = l2(⃗e,⃗e)/√N + 1.
Thel1 orl2 errornormsaregivenbyEq.(4)⃗u=⃗v=⃗eform=1orm=2
respectively. [7]
Question 4. [20 marks] Stellar equilibrium.
Consider the Lane-Emden equation, in the form of the initial value problem
h 0 0 ( x ) + x2 h 0 ( x ) + h ( x ) = 0 h(0) = 1, h0(0) = 0
This equation appears singular at x = 0, but one can use de l’Hoˆpital’s rule or a Taylor expansion of h(x) about x = 0 to show that h00(0) = −1/3. Setting z(x) = h0(x), the above 2nd-order differential equation can be reduced to a system of two 1st-order differential equations for h(x) and z(x):
z0(x)= −13 x=0 −x2z(x)−h(x) x>0
h0(x) = z(x) h(0)=1, z(0)=0
(a) Solve the above 1st-order system numerically, with a 4th-order Runge-Kutta method,
using N + 1 = 101 equidistant points in x ∈ [0, π]. Your code should use three
instances of a vector
tabulate them in your report. [10]
(b) Compute the difference e(x) = hnumerical(x) − hexact(x) between your numerical solution hnumerical(x) and the exact solution hexact(x) = x1 sin x. Output the error values
e(x0), e(x10), e(x20), …, e(x100) to the screen and tabulate them in your report. [5]
3
(c) Compute the error norms:
XN i=0
XN l2(⃗e,⃗e) = |ei|2
|ei|,
implementation of Eq. (4) from Question 2b to compute the norms.] [5]
Question 5. [30 marks] Numerical integration. We wish to compute the definite integral
I=Z1 1 dx −1 1+25×2
numerically and compare to the exact result, Iexact = 25 tan−1(5). (a) Use the composite trapezium rule
Zb XN Δx/2, i=0ori=N b−a f(x)dx≃ wifi, wi= Δx 1≤i≤N−1, Δx= N ,
a i=0
to compute the integral I , using N + 1 = 64 equidistant points in x ∈ [−1, 1]. Use three instances of a vector
values fi = f(xi) and weights wi. [Hint: you may use the function from Question 2a to compute the dot product of the vectors wi and fi.] Output to the screen (and list in
your report) your numerical result Itrapezium and the difference Itrapezium − Iexact. [5]
l1(⃗e,⃗e) =
where ei = e(xi) is the error at each grid-point. [Hint: you may use your C++
i=0
(b) Use the composite Hermite integration rule
Zb XN Δx20 0 f(x)dx≃ wifi+ 12[f(a)−f(b)]
a i=0
with the derivatives f0(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 = 64 equidistant points in x ∈ [−1, 1]. Output to the screen (and list in your report) your numerical result IHermite and the difference IHermite − Iexact.
[5]
onagridofN+1=64pointsxi =−cosθi,whereθi =iπ/N,i=0,1,…,Nto
compute the integral I. [Hint: First compute the values of θi, xi, fi = f(xi) and wi and
store them as vectors. Then, you may use the function from Question 2a to compute
the dot product of the vectors wi and fi.] Output to the screen (and list in your report) your numerical result IClenshawCurtis and the difference IClenshawCurtis − Iexact. [10]
(c) Use the Clenshaw-Curtis quadrature rule
Z b XN 1 , i = 0 or i = N
f(x)dx ≃ wifi, wi = N2 P(N−1)/2 2 cos(2kθ ) , 21− i 1≤i≤N−1
a i=0 N k=1 4k2−1
(d) Compute the integral I using a hit and miss 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.
[10]
4
Question 6. [15 marks] Harmonic oscillator. Consider the first-order system:
dq = p dt
dp = −q dt
(a) Use a 2nd-order Runge-Kutta (RK2) midpoint method to evolve the system, with initial conditions q(0) = 0, p(0) = √2 and time-step Δt = 0.1. Stop the evolution at time
t = 100. Describe how your code works. [8]
(b) Output to the screen (and tabulate in your report) the values of the position q(t),
momentum p(t), energy E(t) = 12 (p2 + q2) and the difference e(t) = E(t) − E(0) for
t = 0, t = 1, t = 10 and t = 100. Is E(t) constant numerically? [7]
End of Paper.
5