联系方式

  • QQ:99515681
  • 邮箱:99515681@qq.com
  • 工作时间:8:00-23:00
  • 微信:codehelp

您当前位置:首页 >> R语言程序R语言程序

日期:2022-04-07 07:54

Main Examination period 2021/22
MTH6150: Numerical Computing in C and C++
Assignment date: Monday 14/3/2022
Submission deadline: Monday 2/5/2022 at 23:59 BST
The coursework is due by Monday, 2nd May 2022 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 your code. You must 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), and include the code and its output
in your report. The C++ code for each question must be pasted into the pdf file of
your report, so that it can be commented on when marking. The code used to answer each
question should also be submitted separately, together with the report, in a single cpp file or
multiple cpp files. You may organise the code in different directories, one for each question.
Late submissions will be treated in accordance with current College regulations. QMPlus
automatically screens for plagiarism. Plagiarism is an assessment offence and carries
severe penalties. In the questions below, use long double if your compiler supports it. If
this is not supported, e.g. if you are using Visual Studio, you may use double.
Please read the instructions above carefully to avoid common mistakes. If in doubt, please
ask. Reports that do not contain C++ code in it (but have code included separately in .cpp
files) are subject to a 25% penalty. Reports that consist only of code and no explanation are
subject to a 40% penalty. C++ code submitted without a report is subject to a 50% penalty.
Reports not accompanied by any C++ code cannot be evaluated and will receive 0 marks.
Only material submitted through QMPlus will be accepted.
1
Coursework [100 marks]
Question 1. [20 marks] Self-consistent iteration.
Solve the transcendental equation
x = e?x
using fixed-point iteration. That is, using the initial guess x0 = 1, obtain a sequence of real
numbers
x1 = e?x0
x2 = e?x1
.
.
.
x50 = e?x49
.
.
.
which tends to the value x∞ = 0.567143 . . ., that is, a root of the function f(x) = x ? e?x.
Formally, the iteration can be written as
xi+1 = e?xi for i = 0, 1, 2, . . . with x0 = 1.
The limit x∞ satisfies f(x∞) = 0.
(a) Write a for loop that performs the above iteration, starting from the initial condition
x0 = 1. Use an if and a break statement to stop the loop when the absolute value of
the difference |xi+1 ? xi| between two consecutive iterations is less than a prescribed
tolerance  = 10?15. Use a long double to store the values of x, and the change in
x, between iterations. Use setprecision(18) to print out the final value of x to 18
digits of accuracy. [10]
(b) In how many iterations did your loop converge? [5]
(c) What is the final error in the above transcendental equation? [Hint: use the final value of
x to compute and display the difference x ? e?x.] Is the error what you expected (i.e. is
it smaller than )? [5]
2
Question 2. [20 marks] Inner products, sums and norms.
(a) Write a function named dot that takes as input two vectors A~ = {A1, ..., An} ∈ Rn and
B~ = {B1, ..., Bn} ∈ Rn (of type valarray) and returns their inner
product
A~ ? B~ = Xn
i=1
AiBi (1)
as a long double number. [5]
(b) Use this function to compute the sum Pn
i=1
1
i2 for n = 106. To do so, you may define a
Euclidean n-vector A={1,1/2,1/3,...,1/1000000} as a valarraydouble> and compute the inner product A~ ? A~. Print the difference A~ ? A~ ? π2/6
between your result and the exact value of the infinite sum on the screen.
Remark: Define π = 3.1415926535897932385 to long double accuracy. [5]
(c) Repeat Question 2(a) using compensated (Kahan) summation and name your function
P
cdot. Use the old dot function and the new cdot function to compute the sum n
i=1 c2 for n = 106 and c = 0.1.
Hint: Define a constant Euclidean n-vector C={c,c,...,c} as a valarraydouble> of size n and compute the inner product C~ ? C~ .
Print the difference between your result and the exact value nc2. [5]
(d) 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 valarray A) const {
which returns the weighted norm
lm(A~) = m
vuutXn
i=1
|Ai|
m =
Xn
i=1
|Ai|
m
1/m
(2)
Use this function object to calculate the norm l2(A~) for the n-vector A~ in Question 2b.
Does the quantity l2(A~)2 equal the inner product A~ ? A~ that you obtained above? [Note:
half marks awarded if you use a regular function instead of a function object.] [5]
3
Question 3. [20 marks] Finite differences.
(a) Write a C++ program that uses finite difference methods to numerically evaluate the
second derivative f00(x) of a function f(x) whose values on a fixed grid of points are
specified f(xi) = fi, i = 0, 1, ..., N. Your code should use valarraydouble> containers to store the values of xi, fi and f00
i . 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 f00
i for f00(xi):
f00
0 = fi+2 ? 2fi+1 + fi
Δx2 if i = 0
f00
i = fi+1 ? 2fi + fi?1
Δx2 if i = 1, 2, ..., N ? 1
f00
N = fi ? 2fi?1 + fi?2
Δx2 if i = N
with Δx = 2/N. Demonstrate that your program works by evaluating the second
derivatives of a known function, f(x) = e?16x2
, with N + 1 = 64 points. Compute the
difference between your numerical derivatives and the known analytical ones:
ei = f00
numerical(xi) ? f00
analytical(xi)
at each grid-point. Output the values xi, fi, f00
i , ei of each valarraydouble> on the screen and tabulate (or plot) them in your report to 10 digits of
accuracy. Comment on whether the error is what you expected. [10]
(b) For the same choice of f(x), demonstrate 1st-order convergence, by showing that, as N
increases, the mean error hei decreases proportionally to Δx = 2/N . You may do so by
tabulating the quantity Nhei for different values of N (e.g. N + 1 = 16, 32, 64, 128)
and checking if this quantity is approximately constant. (Alternatively, you may plot
loghei vs. log N and check if the dependence is linear and if the slope is ?1.)
Here, the mean error hei is defined by
hei = 1
N + 1
X
N
i=0
|ei| = 1
N + 1
l1(~e).
Hint: you may use your code from Question 2d to compute the mean error. [10]
4
Question 4. [20 marks] Numerical integration.
We wish to compute the definite integral
I =
Z b
a
p(b ? x)x dx
numerically, with endpoints a = 0 and b = 4, and compare to the exact result, Iexact = 2π. In
Questions 4a, 4b and 4c below, use instances of a valarray to store the
values of the gridpoints xi, function values fi = f(xi) and numerical integration weights wi.
(a) Use the composite trapezium rule
Z b
a
f(x)dx ' X
N
i=0
wifi = ~w ? ~f, ~w = Δx
2 ? {1, 2, 2, ..., 2, 1}
to compute the integral I, using N + 1 = 64 equidistant points in x ∈ [a, b], that is,
xi = a + iΔx, for i = 0, 1, ..., N, where Δx = b?a
N is the grid spacing. Output your
numerical result Itrapezium and the difference Itrapezium ? Iexact. [Hint: Use the function
dot or cdot, created in Question 2a or 2b, to easily evaluate the inner product ~w ? ~f]. [5]
(b) Use the extended Simpson rule
Z b
a
f(x)dx ' X
N
i=0
wifi, ~w = Δx
48 ? {17, 59, 43, 49, 48, 48, ..., 48, 49, 43, 59, 17}
to compute the integral I, using N + 1 = 64 equidistant points in x ∈ [a, b]. Output
your numerical result ISimpson and the difference ISimpson ? Iexact. [5]
(c) Use the Clenshaw-Curtis quadrature rule
Z b
a
f(x)dx ' X
N
i=0
wifi, wi = b ? a
2
?
( 1
N2 , i = 0 or i = N
2
N

1 ? P(N?1)/2
k=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 in
valarray containers, as shown in the lab. Then, you may use the
function from Question 2a or 2c to compute the inner product ~w ? ~f.
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 Mean Value 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 5. [20 marks] Solitons.
The Korteweg–De Vries (KDV) equation is a non-linear partial differential equation
describing solitons. It can be reduced to a first-order system of ordinary differential equations
of the form:
?
??
??
dq
dt = p
dp
dt = ?V 0
(q)
with a cubic potential V (q) = ?

q3 + 1
2 cq2 + Aq), initial conditions q(0) = ?1
2 c, p(0) = 0,
and parameters A = 0, c = 1.
(a) Use a 2nd-order Runge-Kutta (RK2) method to evolve the system from the initial time
tinitial = 0 until the final time tfinal = 10 with N = 1000 constant time steps
ti = tinitial + iΔt of step size Δt = (tfinal ? tinitial)/N, where i = 0, 1, ..., N.
Use valarray containers to store the values ti, q(ti) = qi,
p(ti) = pi, E(ti) = Ei of the time t, position q(t), momentum p(t) and energy
E(t) = 1
2 p2 + V (q) in each time step, for i = 0, 1, ..., N. Describe how your code works.
Output to the screen (and tabulate in your report) the values of q(t), p(t), and E(t) for
t = 0, t = 1, t = 2, ..., t = 10. How well is E(t) conserved numerically? Is this what
you expected? Why? Can you explain why the energy is or is not numerically
conserved? [5]
(b) Compute the difference e(t) = qnumerical(t) ? qexact(t) between your numerical solution
qnumerical(t) and the exact solution qexact(t) = 1
2 c sech2
(
√c
2 t) for all time steps ti. Output
the error values for t = 0, t = 1, t = 2, ..., t = 10 to the screen and tabulate them in
your report. Comment on whether the error is what you expected. [5]
(c) Use the trapezium rule to repeat parts (a) and (b). Comment again on whether energy is
conserved numerically and why.
Hint: the trapezium rule results in an implicit system:
?
????
????
qi+1 ? qi =
Z ti+1
ti
p(t)dt '
Δt
2 (pi+1 + pi)
pi+1 ? pi = ?
Z ti+1
ti
V (q(t))dt '
Δt
2 [V 0
(qi+1) + V 0
(qi)]
for i = 0, 1, ..., N ? 1. To solve this implicit system for each i, one can perform a
self-consistent iteration (similar to Question 1) at each time step, using a nested for
loop (5-10 self-consistent iterations at each time step would suffice in this case). (That is,
for each i, one can use the initial guess qi+1 = qi and pi+1 = pi on the right hand side of
the above system, evaluate the left hand side, use the new values qi+1 and pi+1 and
substitute them into the right hand side, compute the left hand side again, and so forth,
for 5-10 iterations, then move on to the next time step, as demonstrated in the lab.)
[10]
End of Paper.
6

相关文章

版权所有:留学生编程辅导网 2021,All Rights Reserved 联系方式:QQ:99515681 电子信箱:99515681@qq.com
免责声明:本站部分内容从网络整理而来,只供参考!如有版权问题可联系本站删除。