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

Main Examination period 2019/20

MTH6150: Numerical Computing in C and C++

(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?

## 发表评论