MTH739U / MTH739P: Topics in Scientific Computing 2021/22 – Coursework 2

Part A: Coursework [55 marks]

Question 1.  [9 marks]

(a)  Write a program that numerically evaluates the first and second derivatives using         pseudo-spectral methods of a function f(x) whose values on a fixed grid of points are specified f(xi), i = 0, 1, ..., N, using the method described below. Assume that the    points xi are located at the Chebyshev-Gauss-Lobatto nodes

xi  =          +         zi,    zi  = − cos θi,    θi  =        i = 0, 1, . . . , N          (1)

in the interval x ∈ [a, b]. These points are the extrema of the Chebyshev polynomial of the first kind TN (x). Evaluate the derivatives using the approximations

N

f0 (xi) ' X Df(xj),

j=0

N

f00 (xi) ' X Df(xj)                      (2)

j=0

where Dij are the elements of the (N + 1) × (N + 1) derivative matrices, given by ci( − 1)i+j        i  j

ij    = b         −            i = j = 0                                   (3)

     6               i = j = N

where c0  = cN  = 2 and c1, . . . , cN − 1  = 1. The second derivative matrix can be evaluated by D(2)  = (D(1))2 or, with lower round-off error, by

 ( − 1)i+j       z+zizj −2                   i  0, i  N, i  j

The Chebyshev differentiation matrices can be constructed automatically via the Wolfram Language commands:

D1=NDSolve‘FiniteDifferenceDerivative[Derivative[1],X,

"DifferenceOrder"->"Pseudospectral"]@"DifferentiationMatrix"

D2=NDSolve‘FiniteDifferenceDerivative[Derivative[2],X,

"DifferenceOrder"->"Pseudospectral"]@"DifferentiationMatrix"

which respectively return the matrices D(1) , D(2) for a list X of nodes given by Eq. (1).     [3]

(b)  Demonstrate that your program works by evaluating the first and second derivatives of a known function, such as e −x2  or sin(3x). Graph or tabulate the difference between the

numerical and analytical derivatives, ei  = fumerical(xi) − fnalytical(xi) and

(c)  Show that the difference (measured with an error norm such as l1) between your                     numerical derivatives and the known analytical ones approaches zero as e −N. You might        do this by tabulating or graphing eNl1  ∼ constant or ln l1  ∼ constant − N for several             different values of N, where l1  = P |ei| for the first derivative error norm, and                   similarly for the second derivative.                                                                                        [3]

Remark: For stability and accuracy purposes, it is often beneficial to replace each element Dii of the main diagonal with the negative sum of all matrix elements in the same row,

Dii  = − Pji Dji. This ensures that the differentiation matrix returns zero when multiplying a constant vector.

Question 2.  [11 marks] Implicit and explicit integration of ODEs.

Consider the 2D anisotropic harmonic oscillator, described by the first order system:

  dx

   dy

         =  −

where V(x, y) =  (ωxx2 + ωyy2) is the potential. You may assume that ωx  = 1, ωy  = 2.

(a)  Write the system in the vector form d~u/dt = f~(t, ~u), where ~u := {x, y, px, py}. What is         the vector-valued function f~(t, ~u)? Define this function in Wolfram Language (or a                language of your choice).                                                                                                      [1]

(b)  Use a 2nd order Runge-Kutta (RK2) method that increments the vector ~u in each time  step to evolve the system. Use x(0) = 0, y(0) = 0, px(0) = 1, py(0) = 1 as initial          conditions and Δt = 0.1 as the time-step. Stop the evolution at time t = 100. Describe how your code works. Plot the position components x(t), y(t) as functions of time t.     Plot the energy E(t) =  (p + p) + V(x, y) and the difference δE(t) = E(t) − E(0)

as functions of time. Does the RK2 method conserve energy? Why or why not?                 [5] (c)  Repeat (b) using the trapezium rule

~un+1 − ~un  = Ztn+1  f~(t, ~u(t))dt '  (f~n+1 + f~n)                      (5)

to evolve the system, where ~un  = ~u(tn), f~n  = f~(tn, ~un), and Δt = tn+1 − tn is the (constant) time step. Does the trapezium rule conserve energy? Why or why not?

Remark: this method is implicit, because tn+1 appears on both sides of the evolution              equations. To obtain the future values of the momentum and position, one may use a Do loop at each time step (i.e. a nested Do loop). The loop at each time step can use a RK1 or RK2     step as initial guess, and then iterate the above system of equations 10 times, repeatedly          substituting the last value ~un+1 into the function f~n+1  = f~(tn+1, ~un+1) on right-hand side, to  compute the new value of the left-hand side. That is, at each time step, one can perform the    self-consistent iteration

~u = ~un +  [f~(tn+1, ~u) + f~(tn, ~un)]                               (6)

and stop after 10 iterations.                                                                                                             [5]

Question 3.  [15 marks] Explicit and implicit integration of PDEs.

(a)  Consider the hyperboloidal advection equation

∂tu − (1 − x) ∂xu = 0                                          (7)

for the scalar function u(t, x), with initial data u(0, x) = f(x) for some chosen function f(x), in the domain x ∈ [0, 1]. You may assume that f(x) = e −a[b −ln(1 −x)]2 , with

a = 218 and b = −2/3. Show that uanalytical(t, x) = e −a[b+t −ln(1 −x)]2  satisfies Eq. (7).               (You may use computer algebra to show this analytically).                                                     [1]

(b)  Write a program that uses the method of lines to evolve Eq. (7). Use a pseudo-spectral

method for spatial discretization (for a Chebyshev collocation method, the 1st order       differentiation matrix is given in Question 1). Use an explicit Runge-Kutta method (or   equivalent Taylor method) for time integration. Explain how your programme works.     Plot the solution unumerical(t, x) and the error δu(t) = unumerical(t, x) − uanalytical(t, x) as a function ofx at different instants of time t. What happens if the time-step exceeds the

Courant limit?                                                                                                                             [5] (c)  Plot the energy error δE(t) = E(t) − E(0) as a function of time t, where

E = Z0 1 dx (1 − x)(∂xu)2 .                                       (8)

is the energy. Is energy conserved numerically? Why?

Hint: if using a Chebyshev collocation method with the grid-points of Question 1, you    may use the Clenshaw-Curtis quadrature rule from Coursework 1 to calculate the energy integral.

Hint: stop evaluating the energy before the wave reaches the boundaries, otherwise you          will observe energy loss as the pulse exits the domain.                                                         [4]

(d)  Repeat (b) and (c) using an implicit method (such as the trapezium or Hermite rule) for          time integration. What happens if you increase the time-step? Is energy conserved                  numerically? Why?                                                                                                                [5]

Hint: the trapezium rule applied to the discretized advection equation should yield a system of the form:

~un+1 − ~un     =   Z f~(t, ~u(t))dt '  (f~n+1 + f~n)

where the vectors ~un+1 and f~n+1 represent ui(n)+1  = u(tn+1, xi) and f+1

(9)

= f(tn+1, ui(n)+1)

Question 4.  [10 marks]  Generating random numbers.

For the probability distributions detailed below, construct functions to obtain random numbers from these distributions, and test these functions by creating histograms from a sufficient          number of samples and plotting the histograms together with the given distribution functions.

You are allowed to use any uniform number generator. (Use of off-the-shelf generators for a Normal, Rayleigh or χ2 distribution will lead to half marks).

(a) Uniform distribution over the interval [1, 5].

(b) Uniform distribution over the union of the three intervals [0, 1] ∪ [2, 4] ∪ [5, 6].

(c)  Gaussian distribution with a given mean value μ and variance σ2. Use the Box-Muller

[2]

[2]

method or the Marsaglia polar method or another method of your choice, and test your   function by sampling from a Gaussian with μ = 10 and σ2  = 5. [Half marks awarded if

you use an off-the-shelf function like the Mathematica built-in function

NormalDistribution[μ, σ].]                                                                                                       [2] (d)  Rayleigh distribution with probability density function:

x     x2

σ2

for x ≥ 0, where σ is a parameter provided by the user. Test your function by sampling

from the two distributions obtained by setting σ = 1.0 and σ = 2.5. Use the inversion   method discussed in class to compute the random variates. [Half marks awarded if you

use an off-the-shelf function like the Mathematica built-in function

RayleighDistribution[σ].]

(e)  χ2 distribution with probability density function:

for x ≥ 0. Use the inversion method or another method of your choice. [Half marks awarded if you use an off-the-shelf function like the Mathematica built-in function   ChiSquareDistribution[2].]

Question 5.  [10 marks]  Simple sampling.

The Monte-Carlo estimate of an integral is a random variable, whose mean value approaches the real value of the integral as the number of samples tends to infinity. Write a module that

[2]

[2]

implements a (hit-and-miss or mean-value) Monte-Carlo estimate of the integral

I = Z02 x3 + 2x2 dx.                                               (10)

(a)  Run the code for N = 10n , n = 1, 2, 3, ... uniform random numbers in [0, 2], and obtain            the sequence of Monte-Carlo estimates I1, I2, I3 , . . ..                                                              [3]

(b)  Compute the sequence of absolute errors En  = |I − In| and show that the estimate In  ofI

improves when the number N of samples increases.

(c)  Compute and plot the variance σ2 of the Monte-Carlo estimates as a function ofn. How does the variance change when your estimate is based on more samples?

[3]

[4]

Part B: Programming Project and Report [45 marks]

Answer one of the following questions.

Question 6.  [45 marks] Hyperboloidal wave equation.

(a)  Consider the hyperboloidal wave equation

− (1 + x)∂ψ + (1 − 2x2)∂x∂tψ + (1 − x)x2∂ψ − 2x ∂tψ + x(2 − 3x)∂xψ = 0, (11)

for the scalar function ψ(t, x), with initial data ψ(0, x) = f(x), ψ˙ (0, x) = (1 − x) f0 (x) for some chosen function f(x), in the domain x ∈ [0, 1]. You may assume that                f(x) = e −a[b −ln(1 −x)]2 , with a = 218 and b = −2/3. Show that

ψanalytical(t, x) = e −a[b+t −ln(1 −x)]2  satisfies Eq. (11). (You may use computer algebra to    show this analytically). The above wave equation is second order in space and time.       Define an auxiliary variable π = ∂tψ and use it to reduce the wave equation to a system

first order in time, of the form ∂tu =  u, where u =  and  =  is a

matrix of spatial differentiation operators. What are the operators 1 and 2?                    [2]

(b)  Write a program that uses the method of lines to evolve the first order reduced wave

equation. Use a pseudo-spectral method for spatial discretization (for a Chebyshev collocation method, the differentiation matrices are given in Question 1). Use an    explicit Runge-Kutta method (or equivalent Taylor method) for time integration.     Explain how your programme works. Plot the solution unumerical(t, x) and the error

δu(t) = unumerical(t, x) − uanalytical(t, x) as a function of x at different instants of time t.

What happens if the time-step exceeds the Courant limit?                                                    [10] (c)  Plot the energy error δE(t) = E(t) − E(0) as a function of time t, where

E = Z1 dx (1 + x)(∂tψ)2 + x2 (1 − x)(∂xψ)2 .                      (12)

is the energy. Is energy conserved numerically? Why?

Hint: if using a Chebyshev collocation method with the grid-points of Question 1, you   may use the Clenshaw-Curtis quadrature rule from Coursework 1 to calculate the energy integral.

Hint: stop evaluating the energy before the wave reaches the boundaries, otherwise you         will observe energy loss as the pulse exits the domain.                                                         [8]

(d)  Repeat (b) and (c) using an implicit method (such as the trapezium or Hermite rule) for          time integration. What happens if you increase the time-step? Is energy conserved                  numerically? Why?                                                                                                              [10] Hint: the trapezium rule applied to the discretized advection equation should yield a              system of the form:

 ~un+1 − ~un     =

Z f~(t, ~u(t))dt '  (f~n+1 + f~n)                   (13)

7

where the vectors ~un+1 and f~n+1 represent ui(n)+1  = u(tn+1, xi) and f+1  = f(tn+1, ui(n)+1)

(e)  Write a report describing the two different time evolution methods and details of your     implementation. Compute the Courant limit for the method in parts (b) and (d). What      happens if the time step is too large? Reflect on the relative strengths and weaknesses of

these methods. Back up your conclusions with evidence from specific data such as

number of time-steps needed and total run times.                                                                  [15]