# BEST代写-线上编程学术专家

Best代写-最专业靠谱代写IT | CS | 留学生作业 | 编程代写Java | Python |C/C++ | PHP | Matlab | Assignment Project Homework代写

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

### C++代写 | Main Examination MTH6150: Numerical Computing in C and 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
rn+1 = x ln |rn 1| for n = 0, 1, 2, …
(2)
and stops when |rn+1 rn| <. The function rofx should take x as input and may use
r = 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 iterations
did 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 = {u0, u1, …, uN } R
N+1
~
v = {v , v , …, v } R
0
1
N
product
N
X
~u ~v =
uivi
(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
|uivi|
lm(~u,~v) =
(4)
i=0
Use this function object to calculate the norms l1(~u,~v) and l2(~u,~v) for the 3-vectors in
Question 2a. Does the quantity l2(~u,~v) equal the inner product ~u ~v that you obtained
2
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(xi), i = 0, 1, …, N. Your code should use three instances of a vector<double> to
0
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 ꢀnite differencing to
0
compute an approximation for f (x ):
i
3f(x0) + 4f(x1) f(x2)
+ Ox ) for i = 0
x
f0(x0) =
f0(xi) =
f0(xN ) =
2
f(xi+1) f(x
i1) + Ox ) for i = 1, 2, …, N 1
2
x
f(xN2) 4f(x
N1) + 3f(xN ) + Ox ) for i = N
2
x
with Δx = 2/N. Demonstrate2that your program works by evaluating the derivatives of
a known function, f(x) = e , with N + 1 = 16 points. Compute the difference
between your numerical derivatives and the known analytical ones:
x
ei = fn0umerical(xi) fa0nalytical(xi)
at each grid-point. Output the values ei of this vector<double> on the screen and
tabulate (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
hei = N + 1
|ei|.
i=0
Alternatively, you may use the rms (root mean square) error
v
u
u
t
N
X
1
2
1/2
|ei|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 = l1(~e,~e)/(N + 1) or he i = l2(~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
(
h00(x) + h0(x) + h(x) = 0
2
x
h(0) = 1, h0(0) = 0
This 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
z0(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 x0, x10, x20, …, x100 and h(x0), h(x10), h(x20), …, h(x100) to the screen and
[10]
[5]
(
b) Compute the difference e(x) = hnumerical(x) hexact(x) between your numerical solution
1
x
hnumerical(x) and the exact solution hexact(x) = sin x. Output the error values
e(x0), e(x10), e(x20), …, e(x100) to the screen and tabulate them in your report.
3
(
c) Compute the error norms:
v
u
N
N
X
uX
t
2
l1(~e,~e) =
|ei|, l2(~e,~e) =
|ei|
i=0
i=0
where ei = e(xi) 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, Iexact = tan1(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 ‘
wifi, wi =
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 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
.
[5]
(b) Use the composite Hermite integration rule
Z
N
b
Δx2
X
f(x)dx ‘
wifi + 2 [f0(a) f0(b)]
1
a
i=0
with the derivatives f0(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[1, 1]. Output to the screen (and list in your
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 ‘
wifi, wi =
P
2
N
(N1)/2 2 cos(2i)
1
a
k=1
2
i=0
4k 1
on a grid of N + 1 = 64 points xi =cos θ , where θi = iπ/N, i = 0, 1, …, N to
i
compute the integral I. [Hint: First compute the values of θ , x , f = f(x ) and wi and
i
i
i
i
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
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 IMonteCarlo
and the difference IMonteCarlo Iexact
.
4
Question 6. [15 marks] Harmonic oscillator.
Consider the ꢀrst-order system:
dq
dt = p
dp
dt =q
(
a) Use a 2nd-order Runge-Kutta(RK2) midpoint 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]
[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) for
t = 0, t = 1, t = 10 and t = 100. Is E(t) constant numerically?