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