University of Nottingham
School of Mathematical Sciences
MATH4063 Scientific Computation and C++
Submission Date: Monday 8th January 2024, 15:00 (GMT) Assessed Coursework 2
The following questions are to be used for the coursework assessment in the module MATH4063.
A single zip file containing your answers to the questions below and the code you used to obtain these
answers should be submitted electronically via the MATH4063 Moodle page before the deadline at
the top of this page. You should follow the instructions on the accompanying Coursework Submission
template which is also provided on Moodle. Since this work is assessed, your submission must be
entirely your own work (see the University’s policy on Academic Misconduct).
The style and efficiency of your programs is important. A barely-passing solution will include attempts
to write programs which include some of the correct elements of a solution. A borderline distinction
solution will include working code that neatly and efficiently implements the relevant algorithms, and
that shows evidence of testing.
An indication is given of the weighting of each question by means of a figure enclosed by square
brackets, e.g. [12]. All non-integer calculations should be done in double precision.
Background Material
If you have further questions about this background material, please ask for clarification.
Approximating Systems of Ordinary Differential Equations
Cauchy problems, also known as Initial Value Problems (IVPs), consist of finding solutions to a system
of Ordinary Differential Equations (ODEs), given suitable initial conditions. We will be concerned
with the numerical approximation of the solution to the IVP
du(t)
dt
= f(t,u(t)) for t ∈ [t
0
, T] with u(t
0
) = u
0
, (1)
where f is a sufficiently well-behaved function that maps [t
0
, T) × R
d
to R
d
, the initial condition
u
0 ∈ R
d
is a given vector, and the integer d ≥ 1 is the dimension of the problem. We assume that
f satisfies the Lipschitz condition
kf(t, w) − f(t,u)k ≤ λkw − uk for all w,u ∈ R
d
,
where λ > 0 is a real constant independent of w and u. This condition guarantees that the problem
(1) possesses a unique solution.
We seek an approximation to the solution u(t) of (1) at Nt + 1 evenly spaced time points in the
interval [t
0
, T], so we set
t
n = t
0 + n ∆t for 0 < n ≤ Nt where ∆t = (T − t
0
)/Nt
.
The scalar ∆t is referred to as the time-step. We use a superscript n to denote an approximation to
u(t) at the time points {t
n},
u
n ≈ u(t
n
), for 0 ≤ n ≤ Nt
,
and we are interested in the behaviour of the error e
n = u
n−u(t
n
). We expect this error to decrease
as the step size ∆t tends to 0: the sequence of approximations {u
n} will be generated by a numerical
method, which will be said to be convergent if
lim
∆t→0+
Nt max
n=0
ke
n
k = 0 ,
where k · k is a generic norm on R
d
.
Forward Euler Method
The simplest numerical scheme for the solution of first-order ODEs is the forward Euler method:
u
n+1 = u
n + ∆t f(t
n
,u
n
) for 0 ≤ n < Nt
, (2)
with initial condition u
0 = u(t
0
). If f is analytic, it can be shown that the forward Euler method is
convergent and
E(∆t) := Nt max
n=0
ke
n
k = O(∆t).
Since the error behaves as O((∆t)
p
) where p = 1, the forward Euler method is said to be an order
1 method. This method may suffer from numerical instabilities, hence the step size ∆t must be set
to a sufficiently small value during computations.
Trapezoidal Method
Numerical instabilities can be reduced (and sometimes removed completely) by using an implicit
numerical scheme. One such scheme is the trapezoidal method:
u
n+1 = u
n +
1
2
∆t
f(t
n
,u
n
) + f(t
n+1
,u
n+1)
for 0 ≤ n < Nt
, (3)
with initial condition u
0 = u(t
0
). This method is implicit because it involves f(t
n+1
,u
n+1), which
generates a system of equations which must be solved to compute u
n+1
.
Approximating Partial Differential Equations
In this coursework you will use the finite difference method to approximate the solution of a range
of time-dependent partial differential equations (PDEs), of the form
∂u
∂t = Lu for (x, t) ∈ [xmin, xmax] × [t
0
, T] , (4)
with u(x, t0
) = u
0
(x) for x ∈ [xmin, xmax] ,
where u = u(x, t) is a real function of one spatial coordinate x and a time coordinate t, L is a linear
differential operator involving only derivatives with respect to x, xmin < xmax and T > t0
are all real
numbers, and u
0
is a given real function of x. Throughout these exercises, only Dirichlet boundary
conditions will be considered, imposed at x = xmin and/or x = xmax (as appropriate to the PDE
being approximated).
We seek an approximation to the spatial differential operator Lu of (4) at Nx + 1 evenly spaced
points in the interval [xmin, xmax], so we set
xi = xmin + i ∆x for 0 ≤ i ≤ Nx where ∆x = (xmax − xmin)/Nx .
The scalar ∆x is referred to as the space step. At time t the approximate solution to the PDE is
a vector of values u(t) ∈ R
Nx+1, in which ui(t) ≈ u(xi
, t). In this coursework, the error in this
approximation will be measured only at the final time, t = T, by the discrete norm
E(∆x, ∆t) :=
1
Nx + 1
X
Nx
i=0
(u
Nt
i − u(xi
, T))2
!1
2
, (5)
in which we have used the notation u
n
i
to indicate the approximation to u(xi
, tn
). This can be used
to estimate the order of the approximation.
The approach which will be used is known as the method of lines, in which the differential operator
Lu is approximated at each spatial point xi to generate a vector of right-hand side functions f(t,u(t))
for a system of ODEs of the form (1). To illustrate this we consider two standard PDEs.
2
A Parabolic PDE for Diffusion
The one-dimensional diffusion equation is given by
∂u
∂t = D
∂
2u
∂x2
for (x, t) ∈ [xmin, xmax] × [t
0
, T] , (6)
with the initial condition u(x, t0
) = u
0
(x) for x ∈ [xmin, xmax] and Dirichlet boundary conditions
u(xmin, t) = u−(t) and u(xmax, t) = u+(t), where D > 0 is a given real constant and u− and u+ are
given real functions, which may depend on time. One standard finite difference approximation of the
spatial derivative leads to the semi-discretisation
dui(t)
dt
= D
ui+1(t) − 2ui(t) + ui−1(t)
(∆x)
2
=: fi(t,u(t)), (7)
for i = 1, . . . , Nx−1. The application of Dirichlet boundary conditions involves overwriting the values
of u0(t) and uNx
(t) with, respectively, u−(t) and u+(t), at appropriate times so, for the purposes of
implementation, it can be assumed that fi(t,u(t)) = 0 when i = 1, Nx, for t ∈ [t
0
, T]. This fully
defines the vector f(t,u(t)) in (1), which is combined with the chosen time-stepping method.
For the forward Euler method (2), the fully discrete equations for i = 1, . . . , Nx − 1, i.e. the interior
points, are given by
u
n+1
i = u
n
i + ∆t D
u
n
i+1 − 2u
n
i + u
n
i−1
(∆x)
2
, (8)
in which u
n
i ≈ u(xi
, tn
). For the trapezoidal rule (3), the PDE is approximated at the interior points
by the discrete equations
u
n+1
i = u
n
i +
∆t D
2
u
n
i+1 − 2u
n
i + u
n
i−1
(∆x)
2
+
∆t D
2
u
n+1
i+1 − 2u
n+1
i + u
n+1
i−1
(∆x)
2
. (9)
The values of u
0
i
are provided by the initial conditions and, for Dirichlet boundary conditions, the
equations (8) and (9) are replaced by u
n+1
0 = u−(t
n+1) and u
n+1
Nx = u+(t
n+1) for n = 0, . . . , Nt − 1.
A Hyperbolic Equation for Advection
The one-dimensional constant advection equation is given by
∂u
∂t + v
∂u
∂x = 0 for (x, t) ∈ [xmin, xmax] × [t
0
, T] , (10)
with the initial condition u(x, t0
) = u
0
(x) for x ∈ [xmin, xmax] and Dirichlet boundary conditions
u(xmin, t) = u−(t) if v ≥ 0 or u(xmax, t) = u+(t) if v < 0, where v is a given real constant and
u− and u+ are given real functions, which may depend on time. One standard finite difference
approximation of the spatial derivative leads to the semi-discretisation
dui(t)
dt
= −v
ui(t) − ui−1(t)
∆x
=: fi(t,u(t)), (11)
for i = 1, . . . , Nx−1. The application of Dirichlet boundary conditions involves overwriting the values
of u0(t) or uNx
(t) (depending on the sign of v) with, respectively, u−(t) or u+(t), at appropriate
times. As with the diffusion equation, for the purposes of implementation, it can be assumed that
fi(t,u(t)) = 0 when i = 1 (for v ≥ 0) or i = Nx (for v < 0) for t ∈ [t
0
, T]. A set of fully discrete
equations, analogous to (8) and (9) can be derived in exactly the same way as they were for the
diffusion equation.
3
Materials Provided
You should familiarise yourself with the additional code which has been provided in the folder
Templates/ to perform some of the tasks related to this coursework.
• The abstract class ODEInterface encapsulates an interface to an ODE of the form (1), when
the system consists of a single equation, i.e. d = 1.
• The classes Vector and Matrix are slightly modifiied versions of the classes used in Unit 10
on Iterative Linear Solvers.
• The class UniformGrid1D encapsulates the information and methods needed for constructing,
storing and extracting the spatial discretisation points xi (often referred to as the spatial grid)
for a one-dimensional problem.
• The method GaussianElimination implements the Gaussian elimination algorithm (without
pivoting) for solving a system of linear equations. It uses the Vector and Matrix classes. The
implementation provided is written for general matrices.
• The files plotter.py are Python files provided to help create the plots requested. You do not
have to use them: you may prefer to use alternative graphics tools.
Coursework Questions
In Templates/ you will find a set of folders, one for each question. The folders contain a small
amount of code (.hpp, .cpp and .py files) as well as empty files, which you must edit for the
coursework. You can use any software you want to produce the plots requested below.
You must keep the folder structure and all file names as they are in the templates: the
folder Q1 in your submission, for instance, should be self-contained, and should include all the code
necessary to compile and reproduce your results for Question 1. The template folders may also serve
as a checklist for your submission. As part of your submission, you may also add files to the folders
(for example, new classes, output files, plotting routines, etc.). If you do so, then write a brief
README.txt file, containing a short description of each new file. When you attempt Question 2, use
a new folder and put all the files necessary to produce your results in it; if needed, copy some files
from Q1 to Q2, etc.
This coursework requires you to implement finite difference algorithms for approximating initial and
initial-boundary value problems (IVPs and IBVPs) in an object-oriented manner, then use them to
approximate a range of linear, time-dependent, ordinary and partial differential equations in one space
dimension. Your design choices and your ability to implement classes according to the principles of
object orientation will be assessed throughout this coursework.
1. In this question you will use the forward Euler method to approximate the scalar IVP
du
dt
= a u + sin t for t ∈ [0, T] with u(0) = 0 , (12)
where a is a given real constant. The exact solution to this problem is
u(t) = e
at − a sin t − cost
a
2 + 1
for t ∈ [0, T] .
4
(a) Write an abstract class AbstractODESolver which contains the following members:
• Protected variables for initial and final times
double mFinalTime ;
double mInitialTime ;
• A protected pointer for the ODE system under consideration
ODEInterface * mpODESystem ;
• A protected variable for the current state u
n
double mpState ;
• A protected variable for the time-step size ∆t
double mStepSize ;
• A pure virtual public method
virtual void Solve () = 0;
• Any other member that you choose to implement.
[5]
(b) Write a class LinearODE derived from ODEInterface which:
• Overrides the pure virtual method ComputeF in order to evaluate the right-hand side
of (12).
• Overrides the virtual method ComputeAnalyticSolution in order to compute the
exact solution of (12).
[5]
(c) Write a class ForwardEulerSolver, derived from AbstractODESolver, with the following
members:
• A public constructor
ForwardEulerSolver ( ODEInterface & anODESystem ,
const double initialState ,
const double initialTime ,
const double finalTime ,
const double stepSize ,
const std :: string outputFileName =" output . dat ",
const int saveGap = 1 ,
const int printGap = 1) ;
in which initialState provides the value of u(t
0
).
• A public solution method
void Solve () ;
which computes {u
n} using the forward Euler method for a generic first-order scalar
IVP of the form (1), saves selected elements of the sequences {t
n}, {u
n} in a file, and
prints on screen an initial header and selected elements of the sequences {t
n}, {u
n}.
The method should save to file every saveGap iterations and print on screen every
printGap iterations.
• Any other member that you choose to implement.
[15]
5
(d) Write and execute a main Driver.cpp file which:
i. Approximates the IVP (12) for a = −1, T = 10, using the forward Euler method with
∆t = 0.05, and outputs the solution to a file.
Use your output to plot the approximate solution {u
n} for t ∈ [0, 10], and provide the
approximate value obtained for u(10).
ii. Approximates the IVP (12) with a = −1, T = 1 using the forward Euler method with
various values of ∆t of your choice, computes the corresponding errors E(∆t), and
saves the sequences {∆tk}, {E(∆tk)} to a file.
Use your output to plot log E(∆t) as a function of log ∆t. Include in your report the
values of ∆t and E(∆t) that you used to produce the plot and a brief explanation of
why your results demonstrate that E(∆t) = O(∆t).
Your choices for computing these errors and presenting this evidence will be assessed.
[10]
2. (a) Modify the class ForwardEulerSolver to create a new class TrapezoidalSolver, also
derived from AbstractODESolver, which computes {u
n} using the trapezoidal method
for a generic linear, scalar, first-order IVP of the form (1).
For the purposes of implementation, it is useful to consider the linear ODE in the form
du
dt
= a u + g(t),
for which the two approximations can be written
u
n+1 = u
n + ∆t F(t
n
, un
) Forward Euler (13)
1 −
∆t
2
a
u
n+1 = u
n +
∆t
2
F(t
n
, un
) + ∆t
2
g(t
n+1) Trapezoidal (14)
in which the constant a and the functions F and g depend only on the ODE, not the
discretisation.
You should modify the classes ODEInterface and LinearODE to ensure that your code
retains its encapsulation of the ODE system in this special case. You do not need to
redesign the code to enable it to solve more general ODEs.
[5]
(b) Write and execute a main Driver.cpp file which:
i. Approximates the IVP (12) for a = −1, T = 10, using the trapezoidal method with
∆t = 0.05, and outputs the solution to a file.
Use your output to plot the approximate solution {u
n} for t ∈ [0, 10], and provide the
approximate value obtained for u(10).
ii. Approximates the IVP (12) with a = −1, T = 1 using the trapezoidal method with
various values of ∆t of your choice, computes the corresponding errors E(∆t) and saves
the sequences {∆tk}, {E(∆tk)} to a file.
Use your output to plot log E(∆t) as a function of log ∆t and determine the order of
the method, i.e. the value of p for which E(∆t) = O((∆t)
p
). Include in your report
the values of ∆t and E(∆t) that you used to produce the plot and a brief explanation
of how you determined the value of p.
Is the trapezoidal method better or worse than the forward Euler method for approximating the ODE (12)? Provide a brief justification for your answer.
[5]
6
3. This question concerns the approximation of the one-dimensional diffusion equation using the
methods described in the background material. From the discrete forms (8) and (9), it can be
seen that, for this PDE, the approximations can be written as
u
n+1 = u
n + ∆t F(u
n
) Forward Euler (15)
I −
∆t
2
A
u
n+1 = u
n +
∆t
2
F(u
n
) Trapezoidal (16)
in which the matrix A and the vector F depend only on the discrete form of the spatial operator
Lu. The boundary equations are treated differently from the interior equations because the
Dirichlet boundary conditions are used to overwrite the values of u
n+1
0
and u
n+1
Nx
, so the first
and last rows of A and F need to be defined accordingly.
Note: A good first step for this question would be to copy the relevant code from Q1 and Q2
in to a new folder, convert all double variables used to store values of the approximate solution
u and right-hand side F to Vector variables of length 1, and check that your code still gives
the same answers.
(a) Modify the abstract class AbstractODESolver and the derived classes for the methods
ForwardEulerSolver and TrapezoidalSolver so that the state u(t
n
) is stored in an
object of type Vector. For example, the constructor of the class ForwardEulerSolver
will now take the form
ForwardEulerSolver ( ODEInterface & anODESystem ,
const Vector & initialState ,
const double initialTime ,
const double finalTime ,
const double stepSize ,
const std :: string outputFileName =" output . dat ",
const int saveGap = 1 ,
const int printGap = 1) ;
and the Solve method will have to compute, save and print values of u
n+1 ∈ R
Nx+1
.
• Modify your code so that it computes the discrete norm of the error in the approximation
at the end of the simulation, when t = T, as defined by (5).
• Modify the classes ForwardEulerSolver and TrapezoidalSolver so that they approximate the system of ODEs obtained from the semi-discretisation of a PDE of the
form (4), including the application of Dirichlet boundary conditions.
The trapezoidal method requires the solution of a linear system of equations at each
time-step. You should do this using the method GaussianElimination, which has
been provided. When you are confident that your code is working correctly, you should
modify this method so that it takes full advantage of the tridiagonal structure which the
matrix A has in these cases. Include in your report a brief description of the changes
you have made and the reasons for them.
[10]
(b) Write a class Diffusion, derived from ODEInterface, with the following members:
• A method overriding the method ComputeF of ODEInterface
void ComputeF ( const double t , const Vector & u ,
Vector & f ) const ;
which computes and stores in f the Nx + 1 values of F.
7
• A method overriding the method ComputeAnalyticSolution of ODEInterface
void ComputeAnalyticSolution ( const double t ,
Vector & u ) const ;
which computes the vector of exact solution values u(t) at the points xi
.
• A new method
void ApplyDirichlet ( const double t , Vector & u ) ;
which overwrites the boundary values u0(t) and uNmax (t) of the vector u(t) using the
Dirichlet boundary conditions at the appropriate time level.
• A new method
void ComputeMatrix ( Matrix & A ) const ;
which computes the matrix A.
• Any other method that you choose to implement.
You will need to modify the abstract class ODEInterface to ensure that its design is
consistent with that of the class Diffusion.
[10]
(c) A simple exact solution to the one-dimensional diffusion equation (6) on the interval
[xmin, xmax] = [0, 1], with Dirichlet boundary conditions u+(t) = u−(t) = 0, is
u(x, t) = e−Dπ2
t
sin(πx). (17)
Write and execute a main Driver.cpp file which:
i. Approximates the diffusion equation (6) with D = 0.01 on the interval t ∈ [0, 10] with
initial conditions generated from (17) when t
0 = 0, using the forward Euler method
with Nt = 1000 time-steps and Nx = 100 space steps.
Use your output to plot the initial and final approximate solutions, u
0
and u
Nt
, on the
same graph, and provide the value of the error (5) at the end of the simulation, when
t = 10.
ii. Approximates the diffusion equation (6) with D = 0.01 on the interval t ∈ [0, 10] with
initial conditions generated from (17) when t
0 = 0, using the forward Euler method
with Nx = 100, 200, 400, 800, 1600, space steps. You should use Nt = 1000 with
Nx = 100 and then produce two different sets of results:
• As Nx is increased, increase Nt so that Nt ∝ Nx.
• As Nx is increased, increase Nt so that Nt ∝ Nx
2
.
Use your output to plot log E(∆x, ∆t) as a function of log ∆x in both cases. Use
your results to try to determine the order of the method, p. Include in your report
the values of E(∆x, ∆t) that you used to produce the plots, with a brief explanation
of the behaviour of the errors as Nx and Nt are increased and how you used them to
determine a value for p. Compare this value of p with that observed for the forward
Euler method in Q1(d) and explain any difference between them.
iii. Repeats exercises i. and ii. using the trapezoidal method.
In addition to the plots, output and discussion requested in these exercises, include
in your report a brief explanation of any significant differences between the results
obtained for the two time-stepping methods. Which time-stepping method would you
advise a user to choose for this application? You should consider both stability and
accuracy when determining your answer and briefly justify your choice. You might
consider simulations with other values of Nt and Nx for supporting evidence.
[10]
8
4. This question concerns the approximation of the one-dimensional advection equation using the
methods described in the background material. As with the diffusion equation it is possible to
write the approximations in the forms (15) and (16), though the details are different for the
matrix A, the vector F and the Dirichlet boundary conditions.
(a) Modify the class Diffusion to create a new class Advection, also derived from the abstract class ODEInterface, which encapsulates the system of ODEs which is derived from
the approximation of the one-dimensional advection equation given by (11).
[5]
(b) A simple exact solution to the one-dimensional advection equation (10) on the interval
[xmin, xmax] = [0, 4], with Dirichlet boundary condition u(xmin) = 0 (for v > 0), is
u(x, t) =
cos2
(π(x − vt)) if x − vt ∈ [0.5, 1.5]
0 otherwise.
(18)
Write and execute a main Driver.cpp file which:
i. Approximates the advection equation (10) with v = 2 on the interval t ∈ [0, 1] with
initial conditions generated from (18) when t
0 = 0, using both the forward Euler and
trapezoidal methods with Nt = 1000 time-steps and Nx = 100 space steps.
Use your output to plot the initial and final approximations, u
0
and u
Nt
, on the same
graph (one graph for each method), and provide the value of the error (5) at the end
of the simulation, when t = 1.
ii. Repeats the remaining exercises of Q3(c)ii. and Q3(c)iii. using the same values of Nx
and Nt as were used there.
[5]
(c) The advection equation is generally considered to be more difficult to approximate than the
diffusion equation.
Investigate the behaviour of the methods you have already implemented in the case where
v = −2 (for which the boundary condition is u(xmax) = 0). You do not have to carry out
the same simulations as in part (b) but you should include a brief discussion of the stability
and accuracy of the methods, with appropriate supporting evidence.
Modify your code so that the semi-discretisation in Equation (11) is replaced by
dui(t)
dt
= −v
ui+1(t) − ui−1(t)
2∆x
, (19)
and investigate the behaviour of both time-stepping methods when v = 2. Note that you
will need to use Equation (11) instead of Equation (19) when i = Nx. You do not have to
carry out the same simulations as in part (b) but you should include a brief discussion of
the stability and accuracy of the methods, with appropriate supporting evidence.
[5]
5. This question concerns the approximation of the Black-Scholes equation,
∂u
∂t +
1
2
σ
2x
2
∂
2u
∂x2
+ rx
∂u
∂x − ru = 0 for (x, t) ∈ [xmin, xmax] × [t
0
, T] , (20)
with the final condition u(T, x) = max(x − K, 0) for x ∈ [xmin, xmax] and Dirichlet boundary
conditions u(xmin, t) = 0 and u(xmax, t) = x − Ke−r(T −t)
, where t is time, x is stock price, σ
is volatility, r is risk-free interest rate and K is strike price.
9
This PDE is to be approximated using the semi-discretisation given by
dui(t)
dt
= −
1
2
σ
2xi
2 ui+1(t) − 2ui(t) + ui−1(t)
(∆x)
2
− rxi
ui+1(t) − ui−1(t)
2∆x
+ rui(t), (21)
for i = 1, . . . , Nx − 1. Note that the coefficients of the derivatives depend on x, and one of the
Dirichlet boundary conditions is nonzero and time-dependent.
(a) Modify the class Diffusion (or Advection) to create a new class BlackScholes, also
derived from the abstract class ODEInterface, which encapsulates the system of ODEs
which is derived from the semi-discretisation (21).
This equation is solved backwards in time and you will need to work out how to do this
within the framework you have implemented. In your report, you should briefly describe
how you wrote your code so that the forward Euler and trapezoidal methods step backwards
in time instead of forwards in time.
[5]
(b) Write and execute a main Driver.cpp file which approximates the Black-Scholes equation
(20) with K = 100, r = 0.15, σ = 0.05, on the interval (x, t) ∈ [50, 150] × [0, 1] with
the final condition, given below Equation (20), when T = 1. Use both the forward Euler
method and the trapezoidal method with Nt = 10000 time-steps and Nx = 500 space
steps. You do not have to compute the approximation error in this question.
Use your output to plot the initial and final approximate solutions, u
0
and u
Nt
, on the
same graph (one graph for each method). Which time-stepping method would you advise
a user to choose for this application? You should run additional numerical simulations,
with different values of Nt and Nx, to help you to decide, and use them to justify your
choice. You should also include a brief discussion of why it is appropriate to use the centred
difference approximation for the first derivative in (21), using evidence from the numerical
simulations carried out in previous questions to support your argument.
[5]
The output requested in Questions 1d, 2b, 3c, 4b, 4c and 5b should be included in your submission,
along with any other discussion requested, in the format provided by the solution template file.
10