### C++代写 | Main Examination MTH6150: Numerical Computing in C and C++

这个作业是用C++计算坐标、数值、向量等，并绘制成图形

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 ꢀgures. 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 brieꢁy 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 ꢀle, or multiple cpp

ꢀles compressed into a single zip ꢀle. 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

c

_{b}o_{e}_{t}o_{w}rd_{e}i_{e}n_{n}at_{t}e_{h}_{e}r_{t}∈_{w}_{o}[1_{c},_{o}∞_{o}_{r})_{d}o_{i}_{n}r_{a}a_{t}_{e}ra_{s}d_{y}i_{s}a_{t}l_{e}_{m}to_{s}rt_{,}o_{x}is_{=}e c_{x}o_{(}o_{r}r_{)}d_{,}i_{i}n_{s}a_{g}te_{i}_{v}x_{e}_{n}∈_{b}(_{y}−∞, ∞). The transformationx = 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 speciꢀed tolerance

ꢀ

. The iteration can be performed as a for or while loop that computes the sequence

r

_{n}_{+}_{1 }= x − ln |r_{n}− 1| for n = 0, 1, 2, …(2)

and stops when |r

_{n}_{+}_{1}− r_{n}| < ꢀ. The function rofx should take x as input and may user = 0.8x as initial guess to perform the above iteration, in order to ꢀnd r , r , r , …. The

0

1

2

3

value of x should be ꢀxed 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

^{−}. Run your code for a value of x > 2. In how many iterationsdid it converge? What is the ꢀnal value of r? Substitute r into Eq. (1): is it satisꢀed?

12

[10]

1

Question 2. [10 marks] Inner products.

N+1

and

(of type vector<double>) and returns their inner

(

a) Write a function that takes as input two vectors ~u = {u

_{0}, u_{1}, …, u_{N}} ∈ RN+1

~

v = {v , v , …, v } ∈ R

0

1

N

product

N

X

~u ∙ ~v =

u

_{i}v_{i }(3)

i=0

as a double number. Demonstrate that your program works by computing the inner

product of two real Euclidean 3-vectors, u={2,7,2} and v={3,1,4}.

[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<double> u, const

vector<double> v) const {

which returns the weighted norm

v

u

N

uX

t

m

m/2

|u

_{i}v_{i}|l

_{m}(~u,~v) =(4)

i=0

Use this function object to calculate the norms l

_{1}(~u,~v) and l_{2}(~u,~v) for the 3-vectors inQuestion 2a. Does the quantity l

_{2}(~u,~v) equal the inner product ~u ∙ ~v that you obtained2

above?

[5]

Question 3. [15 marks] Finite differences.

a) Write a C++ program that uses ꢀnite difference methods to numerically evaluate the

rst derivative of a function f(x) whose values on a ꢀxed grid of points are speciꢀed

(

ꢀ

f(x

_{i}), i = 0, 1, …, N. Your code should use three instances of a vector<double> to0

store the values of x

_{i}, f(x_{i}) and f (x_{i}). Assume the grid-points are located atx

_{i }= (2i − N)/N on the interval [−1, 1] and use 2nd order ꢀnite differencing to0

compute an approximation for f (x ):

i

−

3f(x

_{0}) + 4f(x_{1}) − f(x_{2})+ O(Δx ) for i = 0

2Δx

f

^{0}(x_{0}) =f

^{0}(x_{i}) =f

^{0}(x_{N }) =2

^{f}

^{(}

^{x}i+1) − f(x

^{i}

^{−}

^{1}

^{) }+ O(Δx ) for i = 1, 2, …, N − 1

2

2Δx

^{f}

^{(}

^{x}N−2) − 4f(x

^{N}

^{−}

^{1}

^{) }

^{+}

^{ 3}

^{f}

^{(}

^{x}

^{N}

^{ )}+ O(Δx ) for i = N

2

2Δx

with Δx = 2/N. Demonstrate

_{2}that your program works by evaluating the derivatives ofa known function, f(x) = e

^{−}, with N + 1 = 16 points. Compute the differencebetween your numerical derivatives and the known analytical ones:

x

e

_{i }= f_{n}^{0}_{u}_{m}_{e}_{r}_{i}_{c}_{a}_{l}(x_{i}) − f_{a}^{0}_{n}_{a}_{l}_{y}_{t}_{i}_{c}_{a}_{l}(x_{i})at each grid-point. Output the values e

_{i}of this vector<double> on the screen andtabulate (or plot) them in your report.

[8]

2

(

b) For the same choice of f(x), demonstrate 2nd-order convergence, by showing that, as N

2 1/2

increases, the mean error hei or the root mean square error he i decrease

2

−2

2 1/2

proportionally to Δx ∝ N . You may do so by tabulating (or log-log plotting) the

2

2

quantity N hei (or N he i ) for different values of N (e.g. N + 1 = 16, 32, 64, 128).

2

2

2 1/2

Is the quantity N hei (or N he i ) approximately constant?

Here, the mean error hei is deꢀned by

N

X

1

^{h}

^{e}

^{i }

^{=}N + 1

|e

_{i}|.i=0

Alternatively, you may use the rms (root mean square) error

v

u

u

t

N

X

1

2

1/2

|e

_{i}|^{2}.he i

=

N + 1

i=0

Hint: Given the vector ~e, you may use your C++ implementation of Eq. (4) from

√

2

1/2

Question 2 to compute hei = l

_{1}(~e,~e)/(N + 1) or he i = l_{2}(~e,~e)/ N + 1.The l or l error norms are given by Eq. (4) ~u = ~v = ~e for m = 1 or m = 2

1

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) + h^{0}(x) + h(x) = 02

x

h(0) = 1, h

^{0}(0) = 0This equation appears singular at x = 0, but one can use de l’Hoˆpital’s rule or a Taylor

00

0

expansion of h(x) about x = 0 to show that h (0) = −1/3. Setting z(x) = h (x), the above

nd-order differential equation can be reduced to a system of two 1st-order differential

equations for h(x) and z(x):

2

(

1

3

2

x

−

−

x = 0

_{z}0

_{(}

_{x}

_{) }

_{= }

z(x) − h(x) x > 0

0

h (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<double> to store the values of x , h(x ), z(x ). Output the

i

i

i

values x

_{0}, x_{1}_{0}, x_{2}_{0}, …, x_{1}_{0}_{0}and h(x_{0}), h(x_{1}_{0}), h(x_{2}_{0}), …, h(x_{1}_{0}_{0}) to the screen andtabulate them in your report.

[10]

[5]

(

b) Compute the difference e(x) = h

_{n}_{u}_{m}_{e}_{r}_{i}_{c}_{a}_{l}(x) − h_{e}_{x}_{a}_{c}_{t}(x) between your numerical solution1

x

h

_{n}_{u}_{m}_{e}_{r}_{i}_{c}_{a}_{l}(x) and the exact solution h_{e}_{x}_{a}_{c}_{t}(x) = sin x. Output the error valuese(x

_{0}), e(x_{1}_{0}), e(x_{2}_{0}), …, e(x_{1}_{0}_{0}) to the screen and tabulate them in your report.3

(

c) Compute the error norms:

v

u

N

N

X

uX

t

2

l

_{1}(~e,~e) =|e

_{i}|, l_{2}(~e,~e) =|e

_{i}|i=0

i=0

where e

_{i}= e(x_{i}) is the error at each grid-point. [Hint: you may use your C++implementation of Eq. (4) from Question 2b to compute the norms.]

[5]

Question 5. [30 marks] Numerical integration.

We wish to compute the deꢀnite integral

Z

1

1

I =

_{2 }dx

+ 25x

_{−}

_{1 }1

numerically and compare to the exact result, I

_{e}_{x}_{a}_{c}_{t}= tan^{−}^{1}(5).2

5

(

a) Use the composite trapezium rule

(

Z

N

b

X

Δx/2, i = 0 or i = N

_{ ,}_{ Δ}_{x}_{ =}b − a_{, }Δx

1 ≤ i ≤ N − 1

f(x)dx ‘

w

_{i}f_{i}, w_{i}=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<double> to store the values of the gridpoints x , function

i

values f

_{i}= f(x_{i}) and weights w_{i}. [Hint: you may use the function from Question 2a tocompute the dot product of the vectors w

_{i}and f_{i}.] Output to the screen (and list inyour report) your numerical result I

_{t}_{r}_{a}_{p}_{e}_{z}_{i}_{u}_{m}and the difference I_{t}_{r}_{a}_{p}_{e}_{z}_{i}_{u}_{m}− I_{e}_{x}_{a}_{c}_{t }.

[5]

(b) Use the composite Hermite integration rule

Z

N

b

_{Δ}

_{x}2

X

f(x)dx ‘

w

_{i}f_{i }+_{ 2}[f^{0}(a) − f^{0}(b)]1

a

i=0

with the derivatives f

^{0}(x) at x = a and x = b evaluated analytically (and the weights w_{i }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 I

and the difference I

− I

.

exact

[5]

Hermite

Hermite

(

c) Use the Clenshaw-Curtis quadrature rule

(

Z

N

1

N

b

X

,

i = 0 or i = N

,

1 ≤ i ≤ N − 1

2

ꢀ

ꢁ

f(x)dx ‘

w

_{i}f_{i}, w_{i}=P

2

N

(N−1)/2 2 cos(2kθ

_{i})1

−

a

k=1

2

i=0

4k −1

on a grid of N + 1 = 64 points x

_{i}= − cos θ , where θ_{i}= iπ/N, i = 0, 1, …, N toi

compute the integral I. [Hint: First compute the values of θ , x , f = f(x ) and w

_{i}andi

i

i

i

store them as vectors. Then, you may use the function from Question 2a to compute

the dot product of the vectors w

_{i}and f_{i}.] Output to the screen (and list in your report)your numerical result I

_{C}_{l}_{e}_{n}_{s}_{h}_{a}_{w}_{C}_{u}_{r}_{t}_{i}_{s}and the difference I_{C}_{l}_{e}_{n}_{s}_{h}_{a}_{w}_{C}_{u}_{r}_{t}_{i}_{s}− I_{e}_{x}_{a}_{c}_{t }d) Compute the integral I using a hit and miss Monte Carlo method with N = 10000

.

[10]

[10]

(

samples. Output to the screen (and list in your report) your numerical result I

_{M}_{o}_{n}_{t}_{e}_{C}_{a}_{r}_{l}_{o }and the difference I

_{M}_{o}_{n}_{t}_{e}_{C}_{a}_{r}_{l}_{o}− I_{e}_{x}_{a}_{c}_{t }.

4

Question 6. [15 marks] Harmonic oscillator.

Consider the ꢀrst-order system:

dq

_{d}

_{t }= p

dp

_{d}

_{t }= −q

(

a) Use a 2nd-order Runge-Kutta

_{√}(RK2) midpoint method to evolve the system, with initialconditions 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]

[7]

(

b) Output to the screen (and tabulate in your report) the values of the position q(t),

momentum p(t), energy E(t) =

_{ 2}(p + q ) and the difference e(t) = E(t) − E(0) fort = 0, t = 1, t = 10 and t = 100. Is E(t) constant numerically?

1

2 2

End of Paper.

5

## 发表评论