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

这个作业是用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  
cboetowrdeiennattehertwo[1c,oor)doinraaterasdyisatlemtosrt,oxis=e cxo(orr)d,iinsagteivxenb(y−∞,). The transformation  
x = r + ln |r1| .  
(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  
tabulate them in your report.  
[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  
your report) your numerical result Itrapezium and the difference Itrapezium Iexact  
.
[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  
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 ‘  
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?  
1
2 2  
End of Paper.  
5
bestdaixie