Main Examination period 2020/21 MTH6150: Numerical Computing in C and C++ Assignment date: Monday 15/3/2021 Submission deadline: Monday 7/6/2021 at 23:59 The coursework is due byMonday, 7th June 2021 at 23:59 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 must 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 comments 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 an online compiler). 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 directories, 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. 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. Coursework [100 marks] 1 Question 1. [10 marks] Self-consistent iteration. Using a pocket calculator, one may notice that by applying the cosine key repeatedly to the value x0 = 1., one obtains a sequence of real numbers x1 = cos x0 = 0.54030230586814 x2 = cos x1 = 0.857553215846393 . . . x20 = cos x19 = 0.739184399771494 . . . which tends to the value x∞ = 0.739085 . . ., which is the point where the functions x and cos x intersect. The iteration can be written as xn+1 = cos xn for n = 0, 1, 2, . . . with x0 = 1. The limit x∞ satisfies the transcendental equation cos x = x. Write a for or while loop that performs the iteration xn+1 = cos xn starting from the initial condition x0 = 1 and stops when the absolute value of the difference |xn+1 − xn| between two consecutive iterations is less than a prescribed tolerance ² = 10−12. 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− cos x.) [10] 2 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 type vector
) and returns their inner product ~u ∙ ~v = N∑ i=0 uivi (1) as a double number. Demonstrate that your program works by computing the inner products ~u ∙~v, ~u ∙ ~u and ~v ∙~v for the Euclidean 3-vectors, u={1,2,3} and v={4,5,6}. [5] (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 u) const { which returns the weighted norm lm(~v) = m √√√√ N∑ i=0 |vi|m (2) Use this function object to calculate the norms l2(~u) and l2(~v) for the 3-vectors in Question 2a. Do the quantities l2(~u)2 and l2(~v)2 equal the inner products ~u ∙ ~u and ~v ∙ ~v that you obtained above? [Note: half marks awarded if you use a regular function instead of a function object.] [5] 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 to store the values of xi, f(xi) and f ′(xi). 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 for f ′(xi): f ′(x0) = −3f(x0) + 4f(x1)− f(x2) 2Δx + O(Δx2) for i = 0 f ′(xi) = f(xi+1)− f(xi−1) 2Δx + O(Δx2) for i = 1, 2, ..., N − 1 f ′(xN) = f(xN−2)− 4f(xN−1) + 3f(xN ) 2Δx + O(Δx2) for i = N with Δx = 2/N . Demonstrate that your program works by evaluating the derivatives of a known function, f(x) = sin 3x, with N + 1 = 16 points. Compute the difference between your numerical derivatives and the known analytical ones: ei = f ′ numerical(xi)− f ′analytical(xi) at each grid-point. Output the values ei of this vector on the screen and tabulate (or plot) them in your report. [8] 3 (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 〈e〉 = 1 N + 1 N∑ i=0 |ei| = 1 N + 1 l1(~e). [7] Question 4. [20 marks] Stellar equilibrium. Consider the Lane-Emden equation, in the form of the initial value problem{ h′′(x) + 2 x h′(x) + h(x) = 0 h(0) = 1, h′(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 h′′(0) = −1/3. Setting z(x) = h′(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): z′(x) = { −1 3 x = 0 − 2 x z(x)− h(x) x > 0 h′(x) = z(x) h(0) = 1, z(0) = 0 (a) Solve the above 1st-order system numerically, with a 3rd-order Runge-Kutta method, using N + 1 = 101 equidistant points in x ∈ [0, π]. Your code should use three instances of a vector to store the values of xi, h(xi), z(xi). Output the values x0, x10, x20, ..., x100 and h(x0), h(x10), h(x20), ..., h(x100) to the screen and 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) = 1x sin x. Output the error values e(x0), e(x10), e(x20), ..., e(x100) to the screen and tabulate them in your report. [5] (c) Compute the error norms: l1(~e) = N∑ i=0 |ei|, l2(~e) = √√√√ N∑ i=0 |ei|2 where ei = e(xi) is the error at each grid-point. [Hint: you may use your C++ implementation of Eq. (2) from Question 2b to compute the norms.] [5] 4 Question 5. [20 marks] Numerical integration. We wish to compute the definite integral I = ∫ 1 0 1 1 + √ x dx numerically and compare to the exact result, Iexact = 2− log 4. (a) Use the composite trapezium rule∫ b a f(x)dx ' N∑ i=0 wifi, wi = { Δx/2, i = 0 or i = N Δx 1 ≤ i ≤ N − 1 , Δx = b− a N , to compute the integral I , using N + 1 = 64 equidistant points in x ∈ [0, 1]. Use three instances of a vector 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 vectors 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 integration rule∫ b a f(x)dx ' N∑ i=0 wifi + Δx2 12 [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 = 64 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∫ b a f(x)dx ' N∑ i=0 wifi, wi = b− a 2 ∗ { 1 N2 , i = 0 or i = N 2 N ( 1−∑(N−1)/2k=1 2 cos(2kθi)4k2−1 ) 1 ≤ i ≤ N − 1 , on a grid of N + 1 = 64 points xi = (1− 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 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. [5] (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. [5] 5 Question 6. [25 marks] Harmonic oscillator. Consider the first-order system: dq dt = p dp dt = − sin q (a) Use a 4th-order Runge-Kutta (RK4) 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) = 1 2 p2 + 1− cos q and the difference e(t) = E(t)− E(0) for t = 0, t = 5, t = 10, ..., t = 100. To what extend is E(t) conserved numerically? Can you explain why the energy is or is not numerically conserved? [7] (c) Construct a template for the RK4 method and use it to repeat (a) and (b). [10] End of Paper. 6 欢迎咨询51作业君