# Finite Difference Methods for the Black-Scholes Equation

## Full text

(1)

### DIVISION OF APPLIED MATHEMATICS

MÄLARDALEN UNIVERSITY SE-721 23 VÄSTERÅS, SWEDEN

(2)

### Division of Applied Mathematics

Bachelor thesis in mathematics / applied mathematics

Date:

2020-05-29

Project name:

Finite Difference Methods for the Black-Scholes Model

Author(s):

Asima Parveen Saleemi

Supervisor(s): Doghonay Arjmand Reviewer: Milica Rancic Examiner: Siyang Wang Comprising: 15 ECTS credits

(3)

Abstract

Financial engineering problems are of great importance in the academic community and Black-Scholes equation is a revolutionary concept in the modern financial theory. Financial instru-ments such as stocks and derivatives can be evaluated using this model. Option evaluation, is extremely important to trade in the stocks. The numerical solutions of the Black-Scholes equation are used to simulate these options.

In this thesis, the explicit and the implicit Euler methods are used for the approximation of Black-scholes partial differential equation and a second order finite difference scheme is used for the spatial derivatives. These temporal and spatial discretizations are used to gain an insight about the stability properties of the explicit and the implicit methods in general. The numerical results show that the explicit methods have some constraints on the stability, whereas, the implicit Euler method is unconditionally stable. It is also demostrated that both the explicit and the implicit Euler methods are only first order convergent in time and this implies too small step-sizes to achieve a good accuracy.

(4)

## Contents

1 Introduction 3

2 Finite Difference Methods 7 2.1 Forward and Backward Euler’s Methods . . . 7 2.2 Accuracy and Stability . . . 8 2.3 Space Discretization . . . 10 3 A Fully Discrete Approximation of the Black-Scholes Equation 13 3.1 Discretization of Black-Scholes Equation using Explicit Euler Method . . . . 14 3.2 Discretization of Black-Scholes Equation using Implicit Euler Method . . . . 16

4 Numerical Results 18

4.1 Convergence Results . . . 18

5 Conclusion 27

6 Appendix 31

6.1 Numerical Scheme for Explicit Euler . . . 31 6.2 Numerical Scheme for Implicit Euler . . . 34 6.3 Bachelor Degree Objectives . . . 37

(5)

## Introduction

Pricing financial derivatives such as pricing options are of great interest because it can be used to minimize losses caused by price fluctuations of the underlying asset. This protection process is called hedging. There are different types of financial derivatives in the market such as options, forwards, futures and swaps. According to Lesmana and Wang there are two types of options, European options and American options, . European options can only be exercised on the expiry date, while American options can be exercised on or before the expiry date, .

Both traders and investors need a proper method for pricing options to determine a suitable price for a call or a put option. In the last decades there have been numerous attempts in introducing different models for pricing options. A giant leap was taken by Black and Scholes in 1973, , where a mathematical model for calculating a reasonable value for an option was proposed. Another important contribution was made by Robert Merton, who developed the formula purposed by Black and Scholes, .

In 1997, Myron Scholes and Robert Merton received Nobel prize because the theory had a significant impact mainly in financial engineering, and other related fields such as Itô’s calculus, . It introduces a hedging portfolio containing the option itself to share the risk between risky and non-risky assets. Then a strategy is selected to eliminate the possibilities and risks of arbitrage by assuming that the price of risky assets and the bond evolves according to diffusions depending on Itô’s calculus. Further, it derives a valuation equation for options similar to that of Black and Scholes. All three of them came to the conclusion that their formulae leading to a partial differential equation could help to determine an equitable value for a call or a put option.

The partial differential equation introduced by Black and Scholes is in the following form, ∂V ∂ t = − 1 2σ 2S2∂2V ∂ S2 − rS ∂V ∂ S. S> 0,t ∈ (0, T ]. (1.1) Here r and σ are positive constant parameters, referred to risk-free interest rate and market volatility, respectively. Both of these are known functions over the life of the option. The option price is denoted by V , and T is the maturity time of an option. The stock price S is represented by a geometric Brownian motion, which implies that S satisfies the following

(6)

stochastic differential equation

dS= µSdt + σ SdW (1.2) with µ > 0 is the drift rate of the stock and W is standard Brownian motion, .

Equation (1.1) is further equiped with suitable final and boundry conditions i.e., V(S, T ) = max(0, K − S), 0 ≤ S < ∞,

V(S,t) ≈ 0, S→ ∞

V(0,t) = K exp(r(T − t)), 0 ≤ t ≤ T V(S,t) = 0, S → ∞ 0 ≤ t ≤ T. Where K is the strike price.

Upon a suitable change of variable in time, equation (1.1) can be turned into initial bound-ary value problem (note that Black-Scholes equation in its original form requires terminal conditions). The resulting equation ( after the change of variables) is classified as a para-bolic partial differentail equation (PDE). In the literature there are various numerical methods for solving PDEs, in particular parabolic PDEs, . As seen in equation (1.1), there are temporal and spatial differential operators in the problem that need to be discretized for a numerical simulation.

In this thesis, we will restrict our focus to solve the following generalised Black and Scholes equations, which is classified as a parabolic PDE with the given conditions that char-acterize the initial and the boundary conditions of the dependent variables.

∂ u ∂ t − a2(x,t) ∂2u ∂ x2− a1(x,t) ∂ u ∂ x− a0(x,t)u = f (x,t) (1.3) where a0(x,t) = −r(x,t), a1(x,t) = r(x,t) − D(x,t) −1 2σ 2(x,t), a2(x,t) = 1 2σ 2(x,t),

where r and σ are positive constant parameters as stated above and D is the divident yield. This equation models a parabolic PDE and a closed form exact solution is not available due to the fact that the coefficients are variable functions of time and space in general. Therefore, an accurate and efficient approximation of equation (1.3) has been the central goal from a computational point of view. There have been substantial amount of work in the literature

(7)

targetting a numerical approximation of the Black and Scholes equation (1.3), and in general parabolic equations, where the main goal has been to develop computationally efficient and higly accurate numerical methods. A numerical approximation requires discretisations in time and space mainly due to the limited computational power available in the computers.

The Black-Scholes equation was numerically solved in [16, 17] where the authors used the cubic B-spline collocation method and Crank- Nicolson time-stepping scheme in space and time direction respectively to get the second-order accuracy in both space and time for the generalized model. In work  the author obtained third-order accuracy in space using the quantic B-spline collocation method and first and second-order accuracy in time using back-ward Euler method and Crank-Nicolson method, respectively. The WENO method was intro-duced in 1994 for solving convection dominated partial differential equations . Later in 2012, modification of the WENO method extended for the spatial discretization of degenerate parabolic equations which may contain discontinuous solution . In work , the authors have achieved second-order accuracy using an exponential time integration scheme along with a central difference scheme for spatial discretization. In  cubic polynomial spline method and implicit Euler method was used to solve the Black-Scholes model, the authors achieved second-order accuracy in space.

The content of this thesis is not only related to the numerical approximation of the solu-tion to the Black-Scholes model, but it is also relevant for other diffusion problems which are modelled by equation (1.3). Many problems in physics and engineering include instantan-euous changes in time and space, which are often modelled by partial differential equations (PDEs). In nature we are often interested in understanding the behaviour of physical systems whose dynamics change over time within a three-dimensional space. PDEs include equations relating the derivatives of functions of several independent variables, typically representing the spatial and temporal variations. Common application areas include , but are not limited to the wave propagation, heat flow, modelling different physical phenomena from the subatomic structure, e.g., quantum mechanics to problems at larger scales, e.g., continuum mechanics, (see  for a detailed exposition of the theory of PDEs). A class of PDEs, which we often face in applications, can be classified based on the physical nature of the problem they model. For example, hyperbolic equations are used to model the wave propogation, while parabolic and elliptic equations are used to model the diffusion phenomenon in different branches.

Although PDEs are very useful in modelling problems from a broad range of fields, a com-mon difficulty is that exact solutions may be unattainable mainly due to: a) complicated nature of the functions and how their (spatial and temporal) derivatives are related, b) geometrical constraints built in the model. Advances in computer technology during the 20th century have accelerated the development of numerical algorithms approximating the solutions of PDEs. In particular, the main motivation has been to find sufficiently good approximate results with a reasonable amount of computational cost, whenver an exact solution is not in access, see [8, 9, 20] for a few pioneering works in the field.

The goal of this thesis is to explore different ways of approximating Black and Scholes par-tial differenpar-tial equation. In particular, we focus on explicit and implicit Euler discretisations for the temporal discretisations, and use second order finite difference approximations for the spatial derivative. The choice of explicit and implicit Euler for the temporal discretisation is intentional as the motivation is to gain insight about 1) the stability properties of explicit

(8)

and implicit methods in general, 2) the need for developing methods with higher orders of accuracy.

(9)

## Finite Difference Methods

In this thesis, we will focus on finite difference methods for solving partial differential equa-tions, . These techniques include replacing partial derivatives by appropriate difference quotients and the resulting system of equations is solved by standard tools from numerical linear algebra, . Numerical methods are approximation methods, but with good care, the approximations can be made to be so good that the values calculated by an analytical method and by numerical method are not distinguishable.

Finite difference methods are used to obtain approximate solution of differential equations or initial value problems. These methods are based on earlier works of Leonhard Euler from 1768, [4, 10]. It is a remarkable fact that Euler’s idea of looking at discrete approximations came long before the introduction of concept of programming in 1842 by Ada Lovelace, the daughter of Lord Byron.

In the remaining of this section, we will focus on the so-called forward and backward Euler’s approximations for the first order derivatives, and review their properties with respect to accuracy and stability.

### Forward and Backward Euler’s Methods

The basic idea behind Euler’s method is to use the idea of local linearity to join various small lines so that they make up an approximation of the actual curve. Euler’s method is essentially based on the definition of derivative, as we know that dydt denotes the slope of y versus t curve. Once the slope f (t, y(t)) is known we can calculate yt+hby multiplying the slope with equidistant

step size (h). We keep on repeating this method until the time interval (0, T ) is covered. A typical form of the first order ODE is

dy

dt = f (t, y(t)), y(t0) = y0. (2.1) By using Euler’s discretization, we can write

y(t + h) − y(t)

(10)

or

y(t + h) ≈ y(t) + h f (t, y(t)) (2.3) Equation (2.3) is known as the forward Euler’s method or explicit Euler method because it gives a tangent line approximation to the first derivative, that is unknown and written expli-citly in terms of previously calculated values.

One can also notice that Euler’s method is also equivalent to the first two terms of the Taylor’s series approximation, . It is also possible to do higher order discretization for dydt and there are many ways of doing this discretization but all these methods exploit a variant of Taylor’s formula, [12, 23] see also [1, 2, 3], where other novel ideas based on the Taylor’s decomposition but several points are used to develop multi-step numerical methods for ODEs. Another alternative way of approximating (2.1) is

y(t + h) − y(t)

h ≈ f (t + h, y(t + h)) (2.4) i.e.,

y(t + h) ≈ y(t) + h f (t + h, y(t + h)). (2.5) Equation (2.5) is known as Backward Euler’s method or implicit Euler method because both the left and the right hand side of the equation (2.5) include the unknown variable, namely, y(t + h). For nonlinear right hand sides, to find this unknown variable, it is essential to import some other technique, for example the Newton-Raphson method. The implicit Euler’s method is comparatively a stable method due to its nature of finding new values by solving equations that involve not only the current state of time t but also the next state t + h.

### Accuracy and Stability

We consider the local truncation error, i.e., the error which occursin a single step. Local truncation error can be defined as the difference between numerical solution y1and the exact

solution y(t0+ h) at time t0+ h after one step.

The numerical solution can be written as

y1= y0+ h f (t0, y0). (2.6)

Assuming that y(t) is smooth enough, we use the Taylor’s expansion y(t0+ h) = y(t0) + hy0(t0) +

h2 2y

00(t

0) + O(h3). (2.7)

The difference of equations (2.6) and (2.7) can be written as y(t0+ h) − y1=

h2 2y

00(t

(11)

Equation (2.8) illustrates that, for a small step size h, the local truncation error is approximately proportional to h2. The same analysis can be done for the implicit Euler’s method, which again leads to an O(h2) local truncation error. Note that dividing equation (2.7) by h, we can see an O(h) error for approximation of the derivative.

Stability is one of the main concepts regarding numerical methods for solving differential equations. A numerical method is ’stable’ if a small change in the problem causes only a small change in the solution. Furthermore, by stability, we mean how a numerical solution yk grows

as k → ∞. This growth is typically dependent on step size h. In other words, we can say that to obtain reasonable approximations, the step-size in Euler’s method has to be chosen small enough depending on the differential equation.

In this section, we aim at quantifying the condition of the step size. We consider a test equation

y0(t) = −λ y(t) (2.9) with initial condition

y(0) = 1 t∈ [0, T ],

where λ is a constant (complex parameter). The exact solution is y(t) = y0exp (−λ t).

If Re(λ ) > 0, then y(t) → 0 when→ ∞. For explicit Euler method, we have

yk+1= (1 − hλ )yk, (2.10) y1= (1 − hλ )y0, ., ., ., yn= (1 − hλ )ny0.

If Re(λ ) > 0, we must expect that the numerical solution also mimics the behaviour of the exact solution; namely yn→ 0 as n → ∞. This requires that the amplification factor is bounded

by 1, i.e.,

| 1 − hλ |< 1. (2.11) The equation (2.11) can be simplified to h < 2

λ. The explicit Euler method is stable for the test

equation (2.9) if the step size h is smaller than 2

λ.

The motivation to look for implicit methods rather than explicit methods is to construct a more stable algorithm.

We now check the stability of the implicit Euler method using the same test equation (2.9). y0(t) = −λ y(t) (2.12)

(12)

By applying the implicit Euler scheme, we get

yk+1= yk+ (hλ )yk+1, (2.13)

we can solve the test equation for yk+1,

yk+1= 1 1 + hλyk, (2.14) y1= 1 1 + hλy0, (2.15) ., ., ., yn= 1 (1 + hλ )ny0.

If Re(λ ) > 0, then yn→ 0 when n → ∞.

The stability condition requires the amplification factor to be bounded, i.e., 1 1 + hλ < 1 or |1 + hλ | > 1.

This condition is satisfied for any positive step-size h.

Hence, the implicit Euler method is unconditionally stable. The implicit Euler method is computationally more expensive than the explicit Euler method in a single time step, but, it has a better stability property, because it allows us to choose larger step-size. The advantage of choosing larger step-size leads to fewer time steps needed to find the solution y(t) at the final time T .

### Space Discretization

Let us consider the spatial derivatives contained in equation (1.3) α (x)∂

2u

∂ x2+ β (x) ∂ u

∂ x+ γ(x)u = f (x). (2.16) To apply the numerical scheme, we assume boundary conditions as

u(xmin) = uL,

and

(13)

Since we have N points, space is discretized by dividing it in N − 1 intervals, step-size h is taken as

h= xmax− xmin N− 1 .

Now we find the points that lie inbetween xmin and xmax by using the follwing discretization

formula

xi= xmin+ (i − 1)h, i= 1, ...., N. (2.17)

The equation (2.16) stated above includes one second-order derivative which needs to be ap-proximated with the Taylor’s polynomial. The second-order derivative shall be estimated from the sum of forward and backward approximations.

We obtain the following from a forward approximation: f(x + h) = f (x) + h f0(x) +1

2h

2f00(x) + ... (2.18)

For a backward approximation, the following is obtained f(x − h) = f (x) − h f0(x) +1

2h

2f00(x) − ... (2.19)

By combining equations (2.18) and (2.19), we get f(x + h) + f (x − h) = [ f (x) + h f0(x) +1 2h 2f00(x) + ...] + [ f (x) − h f0(x) +1 2h 2f00(x) − ....] =⇒ f (x + h) + f (x − h) = 2 f (x) + h2f00(x) + O(h4). =⇒ f00(x) = f(x + h) − 2 f (x) + f (x − h) h2 − O(h4) h2 . (2.20)

We obtain the following approximation for the second-order derivative as h approaches zero, .

f00(x) = f(x + h) − 2 f (x) + f (x − h) h2 + O(h

2).

(2.21) Let us take the first term in the equation (2.16), which is

α (x)∂

2u

∂ x2 = f (x).

As we have previously calculated the second-order derivative from the sum of forward and backward approximations ,this term can be written as follows

αi  ui+1− 2ui+ ui−1 h2  = f (xi). i= 2, 3...., N − 1. (2.22) For i = 2 α2(u3− 2u2+ u1) = h2f(x2).

(14)

For i = 3 α3(u4− 2u3+ u2) = h2f(x3). . . . For i = N − 1 αN−1(uN+1− 2uN−1+ uN−2) = h2f(xN−1). Where αi= α(xi).

Similarly, if we consider the first-order derivative term in the equation (2.16), we get the following βi  ui+1− ui−1 2h  = f (xi) i= 2, 3...., N − 1, (2.23) where βi= β (xi).

If we insert the equations (2.22) and (2.23) in the equation (2.16), we obtain the following αi  ui+1− 2ui+ ui−1 h2  + βi  ui+1− ui−1 2h  + γiui= f (xi), (2.24) and γi= γ(xi).

(15)

## Black-Scholes Equation

Combining the finite difference approximations (in time and in space) from the previous sec-tion gives the following approximasec-tions of the derivatives.

∂ u ∂ t(xi,tn) = α(xi,tn) ∂2u ∂ x2+ β (xi,tn) ∂ u ∂ x+ γ(xi,tn)u + f (xi,tn). (3.1) u(xi,tn) ∼ ui,n.

The space axis is divided into N number of steps, such that xi= xmin+ (i − 1)h, i= 1, 2....N,

and the time axis is divided into M number of steps, such that tn= (n − 1)∆t, n= 1, 2....M, Forward approximation of the time derivative:

∂ u

∂ t(xi,tn) ≈

ui,n+1− ui,n

∆t . (3.2)

Backward approximation of the time derivative: ∂ u

∂ t(xi,tn) ≈

ui,n− ui,n−1

∆t . (3.3)

Central approximation of the spatial derivative: ∂ u

∂ x(xi,tn) ≈

ui+1,n− ui−1,n

2h . (3.4)

Second derivative (central approximation) of the spatial derivative: ∂2u

∂ x2(xi,tn) ≈

ui+1,n− 2ui,n+ ui−1,n

h2 . (3.5)

These approximations of the derivatives will be used to rewrite the Black-Scholes partial dif-ferential equation (3.1).

(16)

### Euler Method

The Black-Scholes partial differential equation will now be approximated using explicit Euler method. A forward approximation will be used for approximating ∂ u

∂ t and a central

approx-imation will be used for approximating ∂ u

∂ x and the second derivative ∂2u

∂ x2. When inserting the

above stated approximations in the Black-Scholes equation, we get the following ui,n+1− ui,n

∆t = αi,n(

ui+1,n− 2ui,n+ ui−1,n

h2 ) + βi,n(

ui+1,n− ui−1,n

2h ) + γi,nui,n+ fi,n (3.6) This can also be written as

ui,n+1= ui,n+

∆t

h2αi,n(ui+1,n− 2ui,n+ ui−1,n) +

∆t

2hβi,n(ui+1,n− ui−1,n) + ∆tγi,nui,n+ ∆t fi,n. (3.7) i= 2, 3..., N − 1,

and

n= 1, 2, .., M − 1. (3.8) This can be formulated in a matrix form,which will be further used for numerical calculations

un= [u2,n, u3,n, u4,n, ...uN−1,n]T un+1− un ∆t = ( ˜An+ ˜Cn+ ˜Dn)un+ ˜fi,n where ˜ A= 1 h2αi,nA

and we denote matrix A as follows

A=               −2 1 0 0 0 . . . 0 1 −2 1 0 0 . . . 0 0 1 −2 1 0 . . . 0 0 0 1 −2 1 . . . 0 0 0 0 1 −2 1 . . 0 . . . 0 . . . 0 . . . 0 0 0 0 0 0 0 0 1 −2               (3.9)

(17)

and the α matrix can be written as αi,n =               α2,n 0 0 0 0 . . . 0 0 α3,n 0 0 0 . . . 0 0 0 α4,n 0 0 . . . 0 0 0 0 α4,n . . . . 0 0 0 0 0 . . . . 0 . . . 0 . . . 0 . . . 0 0 0 0 0 0 0 0 0 αN−1,n               (3.10) ui,n=               u2,n u3,n u4,n u5,n . . . . uN−1,n               (3.11) also ˜ C= 1 2hβi,nC, where C=               0 1 0 0 0 . . . 0 −1 0 1 0 0 . . . 0 0 −1 0 1 0 . . . 0 0 0 −1 0 1 . . . 0 0 0 0 −1 0 1 . . 0 . . . 0 . . . 0 . . . 0 0 0 0 0 0 0 . −1 0               (3.12)

and the β matrix is

βi,n=               β2,n 0 0 0 0 . . . 0 0 β3,n 0 0 0 . . . 0 0 0 β4,n 0 0 . . . 0 0 0 0 β5,n . . . . 0 0 0 0 0 . . . . 0 . . . 0 . . . 0 . . . 0 0 0 0 0 0 0 0 0 βN−1,n               (3.13)

(18)

where For the term u in equation (3.6), the corresponding matrix is an identity matrix of order (N − 2, N − 2).

where

˜

D= γi,nI(N − 2, N − 2)

and the γi,nmatrix can be written as

γi,n=               γ2,n 0 0 0 0 . . . 0 0 γ3,n 0 0 0 . . . 0 0 0 γ4,n 0 0 . . . 0 0 0 0 γ5,n . . . . 0 0 0 0 0 . . . . 0 . . . 0 . . . 0 . . . 0 0 0 0 0 0 0 0 0 γN−1,n               (3.14) also ˜ fi,n=                f2,n+h12(α2,nuL) − β2,n 2h (u1,n) h2f3,n h2f4,n h2f5,n . . . . fN−1,n+ 1 h2(αN−1,nuR) + βN−1,n 2h (uN,n)                (3.15)

### Euler Method

For the implicit Euler method, a backward approximation will be used to approximate ∂ u ∂ t.

The other derivatives will be approximated the same way as we did for explicit Euler method. By putting the values from the equations (3.2)(3.4) and (3.5) in equation (3.1), we get the following un+1− un ∆t = ( ˜An+1+ ˜Cn+1+ ˜Dn+1)un+1+ ˜fn+1. ⇒ un+1= un+ ∆t( ˜An+1+ ˜Cn+1+ ˜Dn+1)un+1+ ∆t ˜fn+1. [I − ∆t( ˜An+1+ ˜Cn+1+ ˜Dn+1)un+1] = un+ ∆t ˜fn+1. Assume that I− ∆t( ˜An+1+ ˜Cn+1+ ˜Dn+1)un+1= Ln+1,

(19)

where

Ln+1= un+ ∆t ˜fn+1

with the initial condition

u(x, 0) = φ (x); x ∈ [xmin, xmax]

and the boundry conditions

g1(t) ≡ u(xmin,t) = 0;t ∈ [0, T ]

g2(t) ≡ u(xmax,t) = 0;t ∈ [0, T ]

where g1and g2are left and right boundary conditions.

Note that in the implicit Euler method, the inverse L−1n+1needs to be computed at every step. The stability condition for the full time and space discretization using the forward Euler method can be obtained by Von-Neuman analysis, . To see this we consider the case where α is a positive constant and β =γ= f =0 and set the equation (1.3) over the infinite do-main (−∞, ∞) in space. Then we start with the ansatz

uj,n= exp(iωxj) ˆun(ω). (3.16)

from (3.7)

uj,n+1= uj,n+ α

∆t

h2(uj+1,n− 2uj,n+ uj−1,n). (3.17)

Upon plugging (3.16) into (3.17) and letting c = α∆t

h2, we obtain ˆ un+1(ω) = ˆun(ω) (1 + c(2cos(ωh)) − 2) . (3.18) By assuming (1 + c(2cos(ωh)) − 2) = G(ω), (3.19) we obtain ˆ un+1(ω) = G(ω) ˆu(ω). (3.20)

Consequently, we must ensure that the magnitude G(ω) is bounded by 1, i.e., |G(ω)| ≤ 1. This implies that

−1 ≤ (1 + c(2cos(ωh)) − 2) ≤ 1. By simplifying this expression, we get c ≤ 12 or

α∆t h2 ≤ 1 2. ⇒ ∆t ≤ h 2 2α.

Following the same procedure, one can find a similar stability condition when α, β , γ, and f are nonzero. Hence, in general, the stability condition for explicit Euler method is

(20)

## Numerical Results

Figures 4.1 and 4.2 show two simulations of the Black-Scholes model using forward Euler method. We checked two different step-sizes in time (dt) to check the stability properties of explicit Euler method. For a stable solution we have the condition that

dt ≈ C(h2), also dt= T M− 1, and dx= xmax− xmin N . where M ≈C(xmax−xmin)T 2N2.

In Figure 4.1, the parameters are chosen as xmin = −1, xmax = 14 , h = 0.0312 and dt =

0.000146. Note that in the Figures we use the notation dx instead of h. The co-efficients in the model are chosen as

a2(x,t) = t(2 + cos(2πx)) (4.1)

a1(x,t) = t(1 + cos(2πx)) (4.2) a0(x,t) = t(2 + cos(4πx)) (4.3) In simulations, we found that the value of constant at which the solution is stable is 0.15. Figure 4.2 shows that when C is sufficiently small (in our case C = 0.15), we get a stable numerical solution but if C gets larger, the numerical solution becomes unstable.

Figure 4.3 and 4.4 show that implicit Euler method is unconditionally stable.

### Convergence Results

Convergence means that the numerical solution approaches the exact solution on the whole time interval [0, T ] as the mesh size in space and time decreases. It is an essential requirement

(21)

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 -1.5 -1 -0.5 0 0.5 1 1.5

Figure 4.1: Stability of explicit Euler method when N = 40, C = 0.15, dt=0.000146, and dx= 0.0312

.

that a numerical method must satisfy. It implies that by taking sufficiently small step-sizes, we get an accurate solution.

Our theoretical expectation for the decay of error concerning dt and h is maxi,n|uexact(xi,tn) − ui,n| ≤ C(dt + h2)

This demonstrates that the numerical method is first-order accurate in time and second-order accurate in space.

Now we simulate to verify these convergence rates of explicit Euler method both in space and in time. Figure 4.5 and 4.6, show the second and the first order convergence in dx and dt, respectively. Figure 4.7 and 4.8 show the order of dx2for implicit Euler method in space and the order of dt in time.

(22)

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 -6 -4 -2 0 2 4 6 10 135

Figure 4.2: Stability of explicit Euler method when N = 40, C = 0.2, dt = 0.00195 and dx is 0.0312

(23)

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 h -1.5 -1 -0.5 0 0.5 1 1.5 Error Exact Numerical

Figure 4.3: Stability of implicit Euler method when N = 40, C = 0.15, dt = 0.00195 and dx is 0.0312

(24)

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 h -1.5 -1 -0.5 0 0.5 1 1.5 Error Exact Numerical

Figure 4.4: Stability of implicit Euler method when N = 40, C = 0.2, dt = 0.00195 and dx is 0.0312

(25)

10-2 10-1 10-5 10-4 10-3 10-2 10-1

(26)

10-5 10-4 10-3 10-6 10-5 10-4 10-3 10-2 10-1 100

(27)

10-2 10-1 10-5 10-4 10-3 10-2 10-1

(28)

10-3 10-2 10-3

10-2

(29)

## Conclusion

In this thesis, the simulation results gave us an insight about the properties of explicit and im-plicit Euler method. It was demonstrated that exim-plicit Euler method has some constraints on the time step, which is why this method is numerically unstable. On the other hand, implicit Euler method was obsereved to be unconditionally stable. In adition, both explicit and implicit Euler methods are first order convergent in time which particularly means that a good approximation is achieved with a very small step-size. This requires a large number of time-discretization, hence, leading to a larger computation time. This motivates the need for developing higher order numerical methods for Black-Scholes equation.

(30)

## Bibliography

 Doghonay Arjmand. Highly Accurate Difference Schemes for the Numerical Solution of Third-Order Ordinary and Partial Differential Equations. Skolan för datavetenskap och kommunikation, Kungliga Tekniska högskolan, 2010.

 Allaberen Ashyralyev and Doghonay Arjmand. A note on the taylor’s decomposition on four points for a third-order differential equation. Applied mathematics and computation, 188(2):1483–1490, 2007.

 Allaberen Ashyralyev, Doghonay Arjmand, and Muhammet Koksal. Taylor’s decompos-ition on four points for solving third-order linear time-varying systems. Journal of the Franklin Institute, 346(7):651–662, 2009.

 Eduard Batschelet. Über die numerische auflösung von randwertproblemen bei ellipt-ischen partiellen differentialgleichungen. Zeitschrift für angewandte Mathematik und Physik ZAMP, 3(3):165–193, 1952.

 Alain Bensoussan. On the theory of option pricing. Acta Applicandae Mathematica, 2(2):139–158, 1984.

 Fischer Black and Myron Scholes. The pricing of options and corporate liabilities. Journal of political economy, 81(3):637–654, 1973.

 Zhongdi Cen, Anbo Le, and Aimin Xu. Exponential time integration and second-order difference scheme for a generalized black-scholes equation. Journal of Applied Math-ematics, 2012, 2012.

 Richard Courant, Kurt Friedrichs, and Hans Lewy. On the partial difference equations of mathematical physics. IBM journal of Research and Development, 11(2):215–234, 1967.

 G Dahlquist. years of numerical instability, part 1, bit 25 (1985), 188-204. MR 86h, 65006, 33.

 Leonhard Euler. Institutionum calculi integralis volumen primum... 1768.

 Lawrence C Evans. Partial differential equations, volume 19. American Mathematical Soc., 2010.

(31)

 Laurene V Fausett, Laurene V Fausett, Laurene V Fausett, and Laurene V Fausett. Ap-plied numerical analysis using MATLAB, volume 1. Prentice hall Upper Saddle River, NJ:, 1999.

 M Hajipour and A Malek. An efficient high order modified weno scheme for nonlinear parabolic equations. International Journal of Applied Mathematics, 24:443–458, 2011.  Jian Huang and Zhongdi Cen. Cubic spline method for a generalized black-scholes

equation. Mathematical Problems in Engineering, 2014, 2014.

 TE Hull, WH Enright, BM Fellen, and AE Sedgwick. Comparing numerical methods for ordinary differential equations. SIAM Journal on Numerical Analysis, 9(4):603–637, 1972.

 Mohan K Kadalbajoo, Lok Pati Tripathi, and Puneet Arora. A robust nonuniform b-spline collocation method for solving the generalized black–scholes equation. IMA Journal of Numerical Analysis, 34(1):252–278, 2014.

 Mohan K Kadalbajoo, Lok Pati Tripathi, and Alpesh Kumar. A cubic b-spline collocation method for a numerical solution of the generalized black–scholes equation. Mathemat-ical and Computer Modelling, 55(3-4):1483–1505, 2012.

 Stig Larsson and Vidar Thomée. Partial differential equations with numerical methods. 2003.

 Stig Larsson and Vidar Thomée. Partial differential equations with numerical methods, volume 45. Springer Science & Business Media, 2008.

 Peter D Lax and Robert D Richtmyer. Survey of the stability of linear finite difference equations. Communications on pure and applied mathematics, 9(2):267–293, 1956.  Donny C Lesmana and Song Wang. An upwind finite difference method for a

nonlin-ear black–scholes equation governing european option valuation under transaction costs. Applied Mathematics and Computation, 219(16):8811–8828, 2013.

 Xu-Dong Liu, Stanley Osher, Tony Chan, et al. Weighted essentially non-oscillatory schemes. Journal of computational physics, 115(1):200–212, 1994.

 John H Mathews, Kurtis D Fink, et al. Numerical methods using MATLAB, volume 4. Pearson prentice hall Upper Saddle River, NJ, 2004.

 Robert C Merton et al. Theory of rational option pricing. Theory of Valuation, pages 229–288, 1973.

 Reza Mohammadi. Quintic b-spline collocation approach for solving generalized black– scholes equation governing option pricing. Computers & Mathematics with Applications, 69(8):777–797, 2015.

(32)

 Paul Wilmott, Susan Howson, Sam Howison, Jeff Dewynne, et al. The mathematics of financial derivatives: a student introduction. Cambridge university press, 1995.

(33)

## Appendix

### Numerical Scheme for Explicit Euler

c l e a r a l l ; c l c ; c l o s e a l l ; x_min = −1; % i n p u t x_min , x_max = 1 / 4 ; % i n p u t x_min , gamma= @( x , t ) t . * ( 1 + c o s ( 4 * p i *x ) ) ; a l p h a = @( x , t ) t . * ( 2 + c o s ( 2 * p i *x ) ) ; b e t a = @( x , t ) t . * ( 1 + c o s ( 2 * p i *x ) ) ; U e x a c t = @( x , t ) t . * s i n ( 2 * p i *x ) ; U e x a c t _ x = @( x , t ) t . * 2 * p i * c o s ( 2 * p i *x ) ; U e x a c t _ t = @( x , t ) s i n ( 2 * p i *x ) ; U e x a c t _ x x = @( x , t ) − t . * 4 * p i ^2* s i n ( 2 * p i *x ) ; g1 = @( t ) U e x a c t ( x_min , t ) ; g2 = @( t ) U e x a c t ( x_max , t ) ; f x = @( x , t ) U e x a c t _ t ( x , t ) −a l p h a ( x , t ) . * Uexact_xx ( x , t )− b e t a ( x , t ) . . . * U e x a c t _ x ( x , t )−gamma ( x , t ) . * U e x a c t ( x , t ) ; T = 1 ; Cons= 0 . 2 ; %ArrayN = [ 4 0 80 160 320 6 4 0 ] ; %F o r C o n v e r g e n c e i n d t u s e t h i s %ArrayN = [ 1 0 20 40 80 1 6 0 ] ; % F o r c o n v e r g e n c e i n dx u s e t h i s %ArrayM = z e r o s ( 1 , l e n g t h ( ArrayN ) ) ; c o u n t e r = 1 ; %f o r N= ArrayN N= 4 0 ;

(34)

M = c e i l ( T / ( Cons * ( x_max−x_min ) ^ 2 ) *N^ 2 ) ; %N= 8 0 ; %ArrayM ( c o u n t e r ) = M; h = ( x_max−x_min ) / ( N− 1 ) ; % S t e p s i z e h x = z e r o s ( 1 ,N ) ; % M i d d l e p o i n t s f o r i = 1 :N x ( i ) = x_min + ( i −1)*h ; end d t = T / ( M− 1 ) ; % S t e p s i z e h t = z e r o s ( 1 ,M) ; % M i d d l e p o i n t s f o r n = 1 :M t ( n ) = ( n −1)* d t ; end u1 = U e x a c t ( x ( 2 : end − 1 ) , 0 ) ’ ; A = z e r o s ( N−2 ,N− 2 ) ; % D e f i n e M a t r i x "A" A( 1 , 1 ) = − 2 ; A ( 1 , 2 ) = 1 ; A( N−2 ,N−2)= −2; A( N−2 ,N− 3 ) = 1 ; f o r j = 2 : N−3 A( j , j − 1 ) = 1 ; A( j , j )= −2; A( j , j + 1 ) = 1 ; end f o r n = 2 :M, t n = ( n −1)* d t + d t ; %a l p h a * u_xx A l p h a M a t r i x = d i a g ( a l p h a ( x ( 2 : end −1 ) , t n ) ) ; A t i l d a = A l p h a M a t r i x *A; %b e t a * u_x B e t a M a t r i x = d i a g ( b e t a ( x ( 2 : end −1 ) , t n ) ) ; C= z e r o s ( N−2 ,N− 2 ) ; C ( 1 , 1 ) = 0 ; C ( 1 , 2 ) = 1 ;

(35)

C ( N−2 ,N− 2 ) = 0 ; C ( N−2 ,N−3)= −1; B e t a M a t r i x = d i a g ( b e t a ( x ( 2 : end −1 ) , t n ) ) ; f o r j = 2 : N−3; C ( j , j +1)= 1 ; C ( j , j ) = 0 ; C ( j , j −1)= −1; end C t i l d a = B e t a M a t r i x *C ; %Gamma*u

GammaMatrix = d i a g ( gamma ( x ( 2 : end −1 ) , t n ) ) ; D= e y e ( N−2 ,N− 2 ) ;

D t i l d a =GammaMatrix *D; %I m p l e m e n t BCs

b = z e r o s ( N− 2 , 1 ) ;

b ( 1 ) = ( d t / h ^2* a l p h a ( x ( 2 ) , t n ) − b e t a ( x ( 2 ) , t n ) * d t / ( 2 * h ) ) * g1 ( t n ) ;

b ( end ) = ( d t / h ^2* a l p h a ( x ( end −1) , t n ) + b e t a ( x ( end −1) , t n ) * d t / ( 2 * h ) ) . . . * g2 ( t n ) ;

L = D − d t / h ^2* A t i l d a − d t * D t i l d a − d t / ( 2 * h ) * C t i l d a ; u2 = L \ ( u1 + d t * f x ( x ( 2 : end −1) , t n ) ’ + b ) ;

u1=u2 ;

E r r o r ( c o u n t e r , n −1) = max ( max ( a b s ( U e x a c t ( x ( 2 : end −1 ) , t n ) ’ − u2 ) ) ) ; end

p l o t ( x ( 2 : end −1 ) , u2 , ’ bo − ’ , x ( 2 : end −1 ) , U e x a c t ( x ( 2 : end −1 ) ,T ) , ’ r * − ’); % E r r o r 2 ( c o u n t e r ) = max ( a b s ( E r r o r ( c o u n t e r , : ) ) ) ;

% c o u n t e r = c o u n t e r + 1 ;

%l o g l o g ( 2 . / ArrayN , E r r o r , ’ o − − ’ , 2 . / ArrayN , ( 2 . / ArrayN ) . ^ 2 , ’ * − − ’ ) ; x l a b e l ( ’ h ’ ) y l a b e l ( ’ E r r o r ’ ) l e g e n d ( ’ E x a c t ’ , ’ N u m e r i c a l ’ ) ; g r i d on t i t l e ( ’ S t a b i l i t y o f I m p l i c i t E u l e r method w i t h r e s p e c t t o c o n s t a n t C ’ . . . , ’ i n t e r p r e t e r ’ , ’ l a t e x ’ )

% l o g l o g ( ( x_max−x_min ) . / ArrayN , E r r o r 2 , ’ o − − ’ ,( x_max−x_min ) . / ArrayN , %(( x_max−x_min ) . / ArrayN ) . ^ 2 , ’ * − − ’ ) ;

(36)

% x l a b e l ( ’ \$dx\$ ’ , ’ i n t e r p r e t e r ’ , ’ l a t e x ’ ) ; % y l a b e l ( ’ E r r o r ’ , ’ i n t e r p r e t e r ’ , ’ l a t e x ’ ) ; % t i t l e ( ’ C o n v e r g e n c e w i t h r e s p e c t t o \$dx\$ ’ , ’ i n t e r p r e t e r ’ , ’ l a t e x ’ ) ; %l e g e n d ( ’ \$ | u ( x _ i , T ) − u_ { i ,M} | _ { \ i n f t y } \$ ’ , ’ \$O ( dx ^ 2 ) \$ ’ , ’ i n t e r p r e t e r ’ . . . % , ’ l a t e x ’ ) ; s e t ( gca , ’ f o n t s i z e ’ , 1 3 ) ; % g r i d on ; % % f i g u r e ; %

% l o g l o g ( T . / ArrayM , E r r o r 2 , ’ o −−’,T . / ArrayM , ( T . / ArrayM ) , ’ * − − ’); % x l a b e l ( ’ \$ d t \$ ’ , ’ i n t e r p r e t e r ’ , ’ l a t e x ’ ) ; % y l a b e l ( ’ E r r o r ’ , ’ i n t e r p r e t e r ’ , ’ l a t e x ’ ) ; % t i t l e ( ’ C o n v e r g e n c e w i t h r e s p e c t t o \$ d t \$ ’ , ’ i n t e r p r e t e r ’ , ’ l a t e x ’ ) ; %l e g e n d ( ’ \$ | u ( x _ i , T ) − u_ { i ,M} | _ { \ i n f t y } \$ ’ , ’ \$O ( d t ) \$ ’ , ’ i n t e r p r e t e r ’ , . . . %’ l a t e x ’ ) ; % s e t ( gca , ’ f o n t s i z e ’ , 1 3 ) ; % g r i d on ;

### Numerical Scheme for Implicit Euler

c l e a r a l l ; c l c ; c l o s e a l l ; x_min = −1; % i n p u t x_min , x_max = 1 / 4 ; % i n p u t x_min , gamma= @( x , t ) t . * ( 1 + c o s ( 4 * p i *x ) ) ; a l p h a = @( x , t ) t . * ( 2 + c o s ( 2 * p i *x ) ) ; b e t a = @( x , t ) t . * ( 1 + c o s ( 2 * p i *x ) ) ; U e x a c t = @( x , t ) t . * s i n ( 2 * p i *x ) ; U e x a c t _ x = @( x , t ) t . * 2 * p i * c o s ( 2 * p i *x ) ; U e x a c t _ t = @( x , t ) s i n ( 2 * p i *x ) ; U e x a c t _ x x = @( x , t ) − t . * 4 * p i ^2* s i n ( 2 * p i *x ) ; g1 = @( t ) U e x a c t ( x_min , t ) ; g2 = @( t ) U e x a c t ( x_max , t ) ; f x = @( x , t ) U e x a c t _ t ( x , t ) −a l p h a ( x , t ) . * Uexact_xx ( x , t )− b e t a ( x , t ) . . . * U e x a c t _ x ( x , t )−gamma ( x , t ) . * U e x a c t ( x , t ) ; T = 1 ; % ArrayN = [ 4 0 80 160 320 6 4 0 ] ; %F o r C o n v e r g e n c e i n d t u s e t h i s ArrayN = [ 1 0 20 40 80 1 6 0 ] ; % F o r c o n v e r g e n c e i n dx u s e t h i s ArrayM = z e r o s ( 1 , l e n g t h ( ArrayN ) ) ;

(37)

c o u n t e r = 1 ; f o r N= ArrayN %N= 8 0 ; M = 20*N;% c e i l ( T / ( Cons * ( x_max−x_min ) ^ 2 ) *N^ 2 ) ; %N= 8 0 ; ArrayM ( c o u n t e r ) = M; h = ( x_max−x_min ) / ( N− 1 ) ; % S t e p s i z e h x = z e r o s ( 1 ,N ) ; % M i d d l e p o i n t s f o r i = 1 :N x ( i ) = x_min + ( i −1)*h ; end d t = T / ( M− 1 ) ; % S t e p s i z e h t = z e r o s ( 1 ,M) ; % M i d d l e p o i n t s f o r n = 1 :M t ( n ) = ( n −1)* d t ; end u1 = U e x a c t ( x ( 2 : end − 1 ) , 0 ) ’ ; A = z e r o s ( N−2 ,N− 2 ) ; % D e f i n e M a t r i x "A" A( 1 , 1 ) = − 2 ; A ( 1 , 2 ) = 1 ; A( N−2 ,N−2)= −2; A( N−2 ,N− 3 ) = 1 ; f o r j = 2 : N−3 A( j , j − 1 ) = 1 ; A( j , j )= −2; A( j , j + 1 ) = 1 ; end f o r n = 2 :M, t n = ( n −1)* d t + d t ; %a l p h a * u_xx A l p h a M a t r i x = d i a g ( a l p h a ( x ( 2 : end −1 ) , t n ) ) ; A t i l d a = A l p h a M a t r i x *A; %b e t a * u_x B e t a M a t r i x = d i a g ( b e t a ( x ( 2 : end −1 ) , t n ) ) ;

(38)

C= z e r o s ( N−2 ,N− 2 ) ; C ( 1 , 1 ) = 0 ; C ( 1 , 2 ) = 1 ; C ( N−2 ,N− 2 ) = 0 ; C ( N−2 ,N−3)= −1; B e t a M a t r i x = d i a g ( b e t a ( x ( 2 : end −1 ) , t n ) ) ; f o r j = 2 : N−3; C ( j , j +1)= 1 ; C ( j , j ) = 0 ; C ( j , j −1)= −1; end C t i l d a = B e t a M a t r i x *C ; %Gamma*u

GammaMatrix = d i a g ( gamma ( x ( 2 : end −1 ) , t n ) ) ; D= e y e ( N−2 ,N− 2 ) ;

D t i l d a =GammaMatrix *D; %I m p l e m e n t BCs

b = z e r o s ( N− 2 , 1 ) ;

b ( 1 ) = ( d t / h ^2* a l p h a ( x ( 2 ) , t n ) − b e t a ( x ( 2 ) , t n ) * d t / ( 2 * h ) ) * g1 ( t n ) ; b ( end ) = ( d t / h ^2* a l p h a ( x ( end −1) , t n ) + b e t a ( x ( end −1) , t n ) . . .

* d t / ( 2 * h ) ) * g2 ( t n ) ;

L = D − d t / h ^2* A t i l d a − d t * D t i l d a − d t / ( 2 * h ) * C t i l d a ; u2 = L \ ( u1 + d t * f x ( x ( 2 : end −1) , t n ) ’ + b ) ;

u1=u2 ;

E r r o r ( c o u n t e r , n −1) = max ( max ( a b s ( U e x a c t ( x ( 2 : end −1 ) , t n ) ’ − u2 ) ) ) ; end

%p l o t ( x ( 2 : end −1 ) , u2 , ’ bo − ’ , x ( 2 : end −1 ) , U e x a c t ( x ( 2 : end −1 ) ,T ) , ’ r * − ’); E r r o r 2 ( c o u n t e r ) = max ( a b s ( E r r o r ( c o u n t e r , : ) ) ) ;

c o u n t e r = c o u n t e r + 1 ; end

l o g l o g ( ( x_max−x_min ) . / ArrayN , E r r o r 2 , ’ o − − ’ ,( x_max−x_min ) . / ArrayN . . . , ( ( x_max−x_min ) . / ArrayN ) . ^ 2 , ’ * − − ’ ) ;

x l a b e l ( ’ \$dx\$ ’ , ’ i n t e r p r e t e r ’ , ’ l a t e x ’ ) ; y l a b e l ( ’ E r r o r ’ , ’ i n t e r p r e t e r ’ , ’ l a t e x ’ ) ;

(39)

l e g e n d ( ’ \$ | u ( x _ i , T ) − u_ { i ,M} | _ { \ i n f t y } \$ ’ , ’ \$O ( dx ^ 2 ) \$ ’ , ’ i n t e r p r e t e r ’ . . . , ’ l a t e x ’ ) ;

s e t ( gca , ’ f o n t s i z e ’ , 1 3 ) ; g r i d on ;

f i g u r e ;

l o g l o g ( T . / ArrayM , E r r o r 2 , ’ o −−’,T . / ArrayM , ( T . / ArrayM ) , ’ * − − ’); x l a b e l ( ’ \$ d t \$ ’ , ’ i n t e r p r e t e r ’ , ’ l a t e x ’ ) ; y l a b e l ( ’ E r r o r ’ , ’ i n t e r p r e t e r ’ , ’ l a t e x ’ ) ; t i t l e ( ’ C o n v e r g e n c e w i t h r e s p e c t t o \$ d t \$ ’ , ’ i n t e r p r e t e r ’ , ’ l a t e x ’ ) ; l e g e n d ( ’ \$ | u ( x _ i , T ) − u_ { i ,M} | _ { \ i n f t y } \$ ’ , ’ \$O ( d t ) \$ ’ , ’ i n t e r p r e t e r ’ . . . , ’ l a t e x ’ ) ; s e t ( gca , ’ f o n t s i z e ’ , 1 3 ) ; g r i d on ;

### Bachelor Degree Objectives

Here, I provide the objectives from the requirements for obtaining a Bachelor degree from a Swedish University, that I fullfilled throughout my thesis.

• Objective 1: The student should demonstrate knowledge and understanding in the major field of study, including knowledge of the field’s scientific basis, knowledge of applicable methods in the field, specialization in some part of the field and orientation in current research questions.

My thesis report demonstrates this objective as I reflect upon some of the most important concepts in Financial mathematics in Chapter 2. My report specializes in finding the numerical method for solving Black-Scholes equation, where we check the stability and convergence conditions of explicit and implicit Euler method. I also provide orientation in current research question in my Introduction section.

• Objective 2: The student demonstrate the ability to search, collect, evaluate and critic-ally interpret relevant information in a problem formation and to criticcritic-ally discuss the phenomena, problem formulation and situations.

During my work on this report, I went through numerous references to look for a proper and suitable information related to my topic. It was a challenging task to look for right information and to critically evaluate that information to use in the corresponding area. • Objective 3: The student should demonstrate the ability to independently identify,

for-mulate and solve problems and to perform tasks within specified time frames.

I fullfilled this objective by solving the tasks given by my supervisor Doghonay Ar-jmand. He gave me a number of challenging problems to evaluate, which were an essential part of Chapter 3, and I managed to figure out those tasks. Also, there was no delay in terms of deadlines.

(40)

• Objective 4: The student should demonstrate the ability to present orally, in writing and discuss information, problems and solutions in dialogue with dialogue with different groups.

I have a number of constructive discussions with my supervisor from the beginning to the end of thesis writing period. I accomplish this objective on writing abilities throg-hout the report especially in Chapter 1, where I summerize adequately my findings and their implications. I also intend to communicate a clear vision of my thesis during the presentation.

• Objective 5: The student should demonstrate ability in the major field of study, make judgments with respect to scientific, societal and ethical aspects.

This objective is achieved throghout the project as we demonstrated our commitment to use citations correctly. The report respects the area of research and presents opportunit-ies for further research, with an advantage of basic level description of the model.

## Figur Figure 4.1: Stability of explicit Euler method when N = 40, C = 0.15, dt=0.000146, and dx = 0.0312 p.21 Figure 4.3: Stability of implicit Euler method when N = 40, C = 0.15, dt = 0.00195 and dx is 0.0312 p.23 Figure 4.4: Stability of implicit Euler method when N = 40, C = 0.2, dt = 0.00195 and dx is 0.0312 p.24 Figure 4.7: Convergence of implicit euler method in space p.27

Updating...