Numerical Computing in C and C--

MTH6150: Numerical Computing in C and C++

Assignment date: Monday 14/3/2022

Submission deadline: Monday 2/5/2022 at 23:59 BST

The coursework is due by Monday, 2nd May 2022 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 your code. You must 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), and include the code and its output in your report. The C++ code for each question must be pasted into the pdf file of your report, so that it can be commented on when marking. The code used to answer each question should also be submitted separately, together with the report, in a single cpp file or multiple cpp files. You may organise the code in different directories, one for each question. Late submissions will be treated in accordance with current College regulations. QMPlus automatically screens for plagiarism. Plagiarism is an assessment offence and carries severe penalties. In the questions below, use long double if your compiler supports it. If this is not supported, e.g. if you are using Visual Studio, you may use double.

Please read the instructions above carefully to avoid common mistakes. If in doubt, please ask. Reports that do not contain C++ code in it (but have code included separately in .cpp files) are subject to a 25% penalty. Reports that consist only of code and no explanation are subject to a 40% penalty. C++ code submitted without a report is subject to a 50% penalty. Reports not accompanied by any C++ code cannot be evaluated and will receive 0 marks. Only material submitted through QMPlus will be accepted.

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 s = 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]

−

x to compute and display the difference x e−x.] Is the error what you expected (i.e. is

it smaller than s)? [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

A˙ ∙ B˙

= AiBi (1)

i=1

as a long double number. [5]

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

{ }

i=1 i2

∙ ∙ −

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˙ π2/6 between your result and the exact value of the infinite sum on the screen.

Remark: Define π = 3.1415926535897932385 to long double accuracy. [5]

cdot. Use the old dot function and the new cdot function to compute the sum

i=1

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˙) =

‚ n

.Σ

i=1

|Ai|m =

. Σ

i=1

|Ai|

Σ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 jj(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 fijj. 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 fijj for f jj(xi):

f0jj fijj fNjj

= 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 = fnjjumerical(xi) − fajjnalytical(xi)

at each grid-point. Output the values xi, fi, fijj, ei of each valarray<long 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

i=0

N + 1

N + 1 1

(e) = 1 Σ |e | = 1 l (˙e).

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

(b − x)x dx

∫ b √

I =

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.