MTH6150代写、C++程序语言代写

日期：2022-04-06 09:52

Main Examination period 2021/22

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.

1

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]

2

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~ = Xn

i=1

AiBi (1)

as a long double number. [5]

(b) Use this function to compute the sum Pn

i=1

1

i2 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~ − π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]

(c) Repeat Question 2(a) using compensated (Kahan) summation and name your function

P

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

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~) = m

vuutXn

i=1

|Ai|

m =

Xn

i=1

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

3

Question 3. [20 marks] Finite differences.

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

second derivative f00(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 f00

i . 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 f00

i for f00(xi):

f00

0 = fi+2 − 2fi+1 + fi

Δx2 if i = 0

f00

i = fi+1 − 2fi + fi−1

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

f00

N = fi − 2fi−1 + fi−2

Δx2 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 = f00

numerical(xi) − f00

analytical(xi)

at each grid-point. Output the values xi, fi, f00

i , 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 hei decreases proportionally to Δx = 2/N . You may do so by

tabulating the quantity Nhei 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

loghei vs. log N and check if the dependence is linear and if the slope is −1.)

Here, the mean error hei is defined by

hei = 1

N + 1

X

N

i=0

|ei| = 1

N + 1

l1(~e).

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

4

Question 4. [20 marks] Numerical integration.

We wish to compute the definite integral

I =

Z b

a

p(b − x)x dx

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

Z b

a

f(x)dx ' X

N

i=0

wifi = ~w ∙ ~f, ~w = Δx

2 ∗ {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 = b−a

N 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

Z b

a

f(x)dx ' X

N

i=0

wifi, ~w = Δx

48 ∗ {17, 59, 43, 49, 48, 48, ..., 48, 49, 43, 59, 17}

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

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

(c) Use the Clenshaw-Curtis quadrature rule

Z b

a

f(x)dx ' X

N

i=0

wifi, wi = b − a

2

∗

( 1

N2 , i = 0 or i = N

2

N



1 − P(N−1)/2

k=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 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. [5]

(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

Question 5. [20 marks] Solitons.

The Korteweg–De Vries (KDV) equation is a non-linear partial differential equation

describing solitons. It can be reduced to a first-order system of ordinary differential equations

of the form:







dq

dt = p

dp

dt = −V 0

(q)

with a cubic potential V (q) = −

q3 + 1

2 cq2 + Aq), initial conditions q(0) = −1

2 c, p(0) = 0,

and parameters A = 0, c = 1.

(a) Use a 2nd-order Runge-Kutta (RK2) method to evolve the system from the initial time

tinitial = 0 until the final time tfinal = 10 with N = 1000 constant time steps

ti = tinitial + iΔt of step size Δt = (tfinal − tinitial)/N, where i = 0, 1, ..., N.

Use valarray<long double> containers to store the values ti, q(ti) = qi,

p(ti) = pi, E(ti) = Ei of the time t, position q(t), momentum p(t) and energy

E(t) = 1

2 p2 + V (q) in each time step, for i = 0, 1, ..., N. Describe how your code works.

Output to the screen (and tabulate in your report) the values of q(t), p(t), and E(t) for

t = 0, t = 1, t = 2, ..., t = 10. How well is E(t) conserved numerically? Is this what

you expected? Why? Can you explain why the energy is or is not numerically

conserved? [5]

(b) Compute the difference e(t) = qnumerical(t) − qexact(t) between your numerical solution

qnumerical(t) and the exact solution qexact(t) = 1

2 c sech2

(

√c

2 t) for all time steps ti. Output

the error values for t = 0, t = 1, t = 2, ..., t = 10 to the screen and tabulate them in

your report. Comment on whether the error is what you expected. [5]

(c) Use the trapezium rule to repeat parts (a) and (b). Comment again on whether energy is

conserved numerically and why.

Hint: the trapezium rule results in an implicit system:







qi+1 − qi =

Z ti+1

ti

p(t)dt '

Δt

2 (pi+1 + pi)

pi+1 − pi = −

Z ti+1

ti

V (q(t))dt '

Δt

2 [V 0

(qi+1) + V 0

(qi)]

for i = 0, 1, ..., N − 1. To solve this implicit system for each i, one can perform a

self-consistent iteration (similar to Question 1) at each time step, using a nested for

loop (5-10 self-consistent iterations at each time step would suffice in this case). (That is,

for each i, one can use the initial guess qi+1 = qi and pi+1 = pi on the right hand side of

the above system, evaluate the left hand side, use the new values qi+1 and pi+1 and

substitute them into the right hand side, compute the left hand side again, and so forth,

for 5-10 iterations, then move on to the next time step, as demonstrated in the lab.)

[10]

End of Paper.

6