# MATH4063代做、C++编程语言代写

- 首页 >> C/C++编程 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

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