MTH6150: Numerical Computing in C and C++ Main Examination period 2022/23
Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit
Main Examination period 2022/23
MTH6150: Numerical Computing in C and C++
Coursework [100 marks]
Question 1. [20 marks] Self-consistent iteration.
Using a pocket calculator, one may notice that by applying the cosine key repeatedly to the value (in radians) one obtains a sequence of real numbers
x1 = cosx0 = 1.
x2 = cosx1 = 0.54030230586814
x3 = cosx2 = 0.54030230586814
x4 = cosx3 = 0.857553215846393
x21 = cosx20 = 0.739184399771494
which tends to the value x∞ = 0.739085 ..., which is the point where the graphs of the function x and cos x intersect. The iteration can be written as
xn+1 = cosxn for n = 0, 1, 2, . . . with x0 = 1.
The limit x∞ satisfies the transcendental equation
cos x = x.
Write a for loop that performs the iteration xn+1 = cosxn starting from the initial condition x0 = 0 and stops when the absolute value of the difference |xn+1 − xn | between two consecutive iterations is less than a prescribed tolerance e = 10 −12 . Use an If/Break statement to exit the loop when the condition is met. 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 − cosx.) [20]
Question 2. [20 marks] Inner products.
(a) Write a function that takes as input two Euclidean vectors ⃗u = {u1 ,u2 , . . . ,uN } ∈ RN and ⃗v = {v1 ,v2 , . . . ,vN } ∈ RN (of type valarray<long double>) and returns their inner product (also known as Hadamard product)
N
⃗u · ⃗v =工 uivi (1)
i=1
as a long double number. Your function may use (u *v) .sum() to compute the dot product of the valarrays u and v. Create a constant valarray<long double> equal to u={0 . 1,0 . 1, . . . ,0 . 1} with N = 106 elements. Demonstrate
Display the difference ⃗u · ⃗u − 104 on the screen. [5]
(b) Repeat Question 2a using Kahan compensated summation to compute the sum. [5]
(c) 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<double> u) const {
} ...
which returns the weighted norm
「
ℓm (⃗v) = m'/ |vi |m
Use this function object to calculate the norm ℓ2 (⃗u) for the vector in Question 2a. Does the quantity ℓ2 (⃗u)2 equal the inner product ⃗u · ⃗u that you obtained above? [Note: half marks awarded if you use a regular function instead of a function object.] [10]
Question 3. [20 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 valarray<long double> to store the values of xi , f(xi ) and f\ (xi ). Assume the grid-points are located at xi = a + i∆x with ∆x = (b − a)/N . on the interval x ∈ [a,b] and use 2nd order finite differencing to compute an approximation for f\ (xi ):
−3f(x0 ) + 4f(x1 ) − f(x2 )
f\ (xi ) = + O(∆x2 ) for i = 1, 2, . . . ,N − 1
f(xN−2) − 4f(xN−1) + 3f(xN )
2∆x
Demonstrate that your program works by evaluating the derivatives of a known function, f(x) = sin3x, with N + 1 = 32 points on the interval x ∈ [a,b] = [−1, 1]. Compute the difference between your numerical derivatives and the known analytical ones:
ei = fumerical (xi ) − fnalytical (xi )
at each grid-point. Output the values ei of this valarray<long double> on the screen and tabulate (or plot) them in your report. [10]
(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
N
⟨e⟩ = |ei | = ℓ 1 (⃗e).
You may use your code from Question 3c to compute the ℓ1 norm. [10]
Question 4. [20 marks] Numerical integration.
We wish to compute the definite integral
I = \a b sin ( ) dx
numerically for a = 0, b = 10 and compare to the exact result, Iexact = 2.74324739415100920.
(a) Use the composite trapezium rule
\a b f(x)dx ≃ wi fi , wi = ∆x = ,
to compute the integral I, using N + 1 = 128 equidistant points in x ∈ [a,b]. Use three instances of a valarray<long double> 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 valarrays 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 rule
\a b f(x)dx ≃ wi fi + [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 = 128 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
\a b f(x)dx ≃ wi fi ,
wi = ∗ { 2 ,(1 −z1(−)1)/2
i = 0 or i = N
) 1 ≤ i ≤ N − 1 ,
on a grid of N + 1 = 128 points xi = [(a + b) − (b − a)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 valarrays. Then, you may use the function from Question 2a to compute the dot product of the valarrays wi and fi .] Output to the screen (and list in your report) your numerical result IClenshawCurtis and the difference
(d) Compute the integral I using a Mean Value Monte Carlo method with N = 1000,
N = 10000 and N = 100000 samples. Output to the screen (and list in your report)
your numerical results IMonteCarlo and the difference IMonteCarlo − Iexact for each N. [5]
Question 5. [20 marks] Stellar equilibrium.
Consider the generalized Lane-Emden equation, in the form of the initial value problem
2x, (3)
This equation describes the structure of a self-gravitating planet or stellar object with mass density ρ(h) = h, where h is the specific enthalpy. Setting m(x) := −x2 h\ (x), the above 2nd-order differential equation can be reduced to a system of two 1st-order differential equations for h(x) and m(x):
h\ (x) =
〈 m\ (x) = h(x) x2
h(0) = 1
(m(0) = 0
The above equation appears singular at x = 0, but one can use a Taylor expansion about the origin to show that m(x) ≃ x3 + O(x5 ) for small x, so that h\ (0) = 0.
Declare a function
valarray<long double> F(const long double t, const
valarray<long double>& u) {
} ...
that takes the radius x and a valarray with elements ⃗u = {h,m} as arguments and returns a valarray with components f⃗= {h\ ,m\ } as output, with h\ ,m\ given by Eq. (4)
above. [5]
(b) Solve the above 1st-order system numerically, with a 4th-order Runge-Kutta method
long double RK4(const double t, const double dt, const valarray<long double> &, long double f(const long double,
valarray<long double> &) {
} ...
using N + 1 = 151 equidistant points in x ∈ [0,π].
Declare a valarray of valarrays:
valarray<valarray<long double>> u
and use it to store ⃗u(xi ) = {hi ,mi } = {h(xi ),m(xi )} at each grid-point xi .
Output (only) the values {x0 ,x10 ,x20 , . . . ,xN } and {h(x0 ),h(x10 ),h(x20 ), . . . ,h(xN )} to the screen and tabulate them in your report. [5]
(c) Compute the difference e(x) = hnumerical (x) − hexact (x) between your numerical solution hnumerical (x) and the exact solution
hexact (x) = sinc x = {
x 0
x = 0
Output the error values e(x0 ),e(x10 ),e(x20 ), . . . ,e(xN ) to the screen and tabulate them in your report. Comment on whether the error is what you expected and why. What is the radius x = R of the star, where hexact (x) = 0? [Hint: use your answer from Question 1].
Remark: it is more convenient to display {xi }, {h(xi )}, {e(xi )} all at once, that is, answer Questions 5b and 5c in the same loop. [5]
(d) Compute the error norms:
‘
N 图 N
l1 (⃗e) = |ei |, l2 (⃗e) = /图 |ei |2
implementation of Eq. (2) from Question 2 to compute the norms.] Comment on
whether the error norms are what you expected. [5]
2023-04-17