程序代写案例-MTH6150

Main Examination period 2020/21
MTH6150: Numerical Computing in C and C++
Assignment date: Monday 15/3/2021
Submission deadline: Monday 7/6/2021 a
t 23:59
The coursework is due byMonday, 7th June 2021 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 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 an online compiler). The code produced to answer
each question should be submitted aside with the report, in a single cpp file, or multiple cpp
files compressed into a single zip file. It is useful to organise the code in different directories,
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. 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.
Coursework [100 marks]
1
Question 1. [10 marks] Self-consistent iteration.
Using a pocket calculator, one may notice that by applying the cosine key repeatedly to the
value x0 = 1., one obtains a sequence of real numbers
x1 = cos x0 = 0.54030230586814
x2 = cos x1 = 0.857553215846393
.
.
.
x20 = cos x19 = 0.739184399771494
.
.
.
which tends to the value x∞ = 0.739085 . . ., which is the point where the functions x and
cos x intersect. The iteration can be written as
xn+1 = cos xn for n = 0, 1, 2, . . . with x0 = 1.
The limit x∞ satisfies the transcendental equation
cos x = x.
Write a for or while loop that performs the iteration xn+1 = cos xn starting from the initial
condition x0 = 1 and stops when the absolute value of the difference |xn+1 − xn| between two
consecutive iterations is less than a prescribed tolerance ² = 10−12. Print out the final value
xn+1 to 16 digits of accuracy. In how many iterations did your loop converge? What is the
final error in the above transcedental equation? (Hint: use the final value to compute and print
out the difference x− cos x.) [10]
2
Question 2. [10 marks] Inner products.
(a) Write a function that takes as input two vectors ~u = {u0, u1, ..., uN} ∈ RN+1 and
~v = {v0, v1, ..., vN} ∈ RN+1 (of type vector) and returns their inner
product
~u ∙ ~v =
N∑
i=0
uivi (1)
as a double number. Demonstrate that your program works by computing the inner
products ~u ∙~v, ~u ∙ ~u and ~v ∙~v for the Euclidean 3-vectors, u={1,2,3} and v={4,5,6}. [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 u) const {
which returns the weighted norm
lm(~v) =
m
√√√√ N∑
i=0
|vi|m (2)
Use this function object to calculate the norms l2(~u) and l2(~v) for the 3-vectors in
Question 2a. Do the quantities l2(~u)2 and l2(~v)2 equal the inner products ~u ∙ ~u and ~v ∙ ~v
that you obtained above? [Note: half marks awarded if you use a regular function
instead of a function object.] [5]
Question 3. [15 marks] Finite differences.
(a) Write a C++ program that uses finite difference methods to numerically evaluate the
first derivative of a function f(x) whose values on a fixed grid of points are specified
f(xi), i = 0, 1, ..., N . Your code should use three instances of a vector to
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 finite differencing to
compute an approximation for f ′(xi):
f ′(x0) =
−3f(x0) + 4f(x1)− f(x2)
2Δx
+ O(Δx2) for i = 0
f ′(xi) =
f(xi+1)− f(xi−1)
2Δx
+ O(Δx2) for i = 1, 2, ..., N − 1
f ′(xN) =
f(xN−2)− 4f(xN−1) + 3f(xN )
2Δx
+ O(Δx2) for i = N
with Δx = 2/N . Demonstrate that your program works by evaluating the derivatives of
a known function, f(x) = sin 3x, with N + 1 = 16 points. Compute the difference
between your numerical derivatives and the known analytical ones:
ei = f

numerical(xi)− f ′analytical(xi)
at each grid-point. Output the values ei of this vector on the screen and
tabulate (or plot) them in your report. [8]
3
(b) For the same choice of f(x), demonstrate 2nd-order convergence, by showing that, as N
increases, the mean error 〈e〉 decreases proportionally to Δx2 ∝ N−2 . You may do so
by tabulating the quantity N2〈e〉 for different values of N (e.g. N + 1 = 16, 32,
64, 128) and checking if this quantity is roughly constant. Alternatively, you may plot
log〈e〉 vs. log N and check if the dependence is linear and if the slope is -2.
Here, the mean error 〈e〉 is defined by
〈e〉 = 1
N + 1
N∑
i=0
|ei| = 1
N + 1
l1(~e).
[7]
Question 4. [20 marks] Stellar equilibrium.
Consider the Lane-Emden equation, in the form of the initial value problem{
h′′(x) + 2
x
h′(x) + h(x) = 0
h(0) = 1, h′(0) = 0
This equation appears singular at x = 0, but one can use de l’Hoˆpital’s rule or a Taylor
expansion of h(x) about x = 0 to show that h′′(0) = −1/3. Setting z(x) = h′(x), the above
2nd-order differential equation can be reduced to a system of two 1st-order differential
equations for h(x) and z(x):
z′(x) =
{
−1
3
x = 0
− 2
x
z(x)− h(x) x > 0
h′(x) = z(x)
h(0) = 1, z(0) = 0
(a) Solve the above 1st-order system numerically, with a 3rd-order Runge-Kutta method,
using N + 1 = 101 equidistant points in x ∈ [0, π]. Your code should use three
instances of a vector to store the values of xi, h(xi), z(xi). Output the
values x0, x10, x20, ..., x100 and h(x0), h(x10), h(x20), ..., h(x100) to the screen and
tabulate them in your report. [10]
(b) Compute the difference e(x) = hnumerical(x)− hexact(x) between your numerical solution
hnumerical(x) and the exact solution hexact(x) = 1x sin x. Output the error values
e(x0), e(x10), e(x20), ..., e(x100) to the screen and tabulate them in your report. [5]
(c) Compute the error norms:
l1(~e) =
N∑
i=0
|ei|, l2(~e) =
√√√√ N∑
i=0
|ei|2
where ei = e(xi) is the error at each grid-point. [Hint: you may use your C++
implementation of Eq. (2) from Question 2b to compute the norms.] [5]
4
Question 5. [20 marks] Numerical integration.
We wish to compute the definite integral
I =
∫ 1
0
1
1 +

x
dx
numerically and compare to the exact result, Iexact = 2− log 4.
(a) Use the composite trapezium rule∫ b
a
f(x)dx '
N∑
i=0
wifi, wi =
{
Δx/2, i = 0 or i = N
Δx 1 ≤ i ≤ N − 1 , Δx =
b− a
N
,
to compute the integral I , using N + 1 = 64 equidistant points in x ∈ [0, 1]. Use three
instances of a vector to store the values of the gridpoints xi, function
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
your report) your numerical result Itrapezium and the difference Itrapezium − Iexact. [5]
(b) Use the composite Hermite integration rule∫ b
a
f(x)dx '
N∑
i=0
wifi +
Δx2
12
[f ′(a)− f ′(b)]
with the derivatives f ′(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 ∈ [0, 1]. Output to the screen (and list in your report)
your numerical result IHermite and the difference IHermite − Iexact. [5]
(c) Use the Clenshaw-Curtis quadrature rule∫ b
a
f(x)dx '
N∑
i=0
wifi, wi =
b− a
2

{
1
N2
, i = 0 or i = N
2
N
(
1−∑(N−1)/2k=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 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. [5]
(d) Compute the integral I using a hit and miss 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 6. [25 marks] Harmonic oscillator.
Consider the first-order system: 
dq
dt
= p
dp
dt
= − sin q
(a) Use a 4th-order Runge-Kutta (RK4) 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]
(b) Output to the screen (and tabulate in your report) the values of the position q(t),
momentum p(t), energy E(t) = 1
2
p2 + 1− cos q and the difference e(t) = E(t)− E(0)
for t = 0, t = 5, t = 10, ..., t = 100. To what extend is E(t) conserved numerically?
Can you explain why the energy is or is not numerically conserved? [7]
(c) Construct a template for the RK4 method and use it to repeat (a) and (b). [10]
End of Paper.
6

欢迎咨询51作业君
51作业君 51作业君

Email:51zuoyejun

@gmail.com

添加客服微信: ITCSdaixie