LGS MF2 Computer Lab Session 1_Explicit

Finite-Difference Schemes for the 
Diffusion Equation with Smooth 
Initial Conditions

Schemes Investigated

In these first sessions we compare the accuracy of various difference schemes for solving the diffusion equation. This is the equation that arises when the Black-Scholes differential equation is transformed into a form suitable for treatment by finite-difference methods. We compare

(a) explicit finite-difference, with 400 time-steps (equivalent to the use of a binomial model, but on a grid rather than a tree);
(b) fully implicit, also with 400 time-steps;
(c) Crank-Nicholson, with 40 time-steps;
(d) Douglas, with 40 time-steps.

The solution method for type (a) is a simple updating rule, while (b), (c), (d) require the solution of tridiagonal systems of equations. This notebook and the associated C++ code look at case (a) above.

A Simple Test Problem with Smooth Initial Conditions

We consider the diffusion equation

((∂u)/(∂τ))=(((∂^2)u)/(∂(x^2)))

on the region defined by

-2≤x≤2          τ≥0

The initial condition is

u(x,0)=sin(((πx)/2))

and the boundary conditions are

u(2,τ)=u(-2,τ)=0

This has the exact solution

u(x,τ)=sin(((πx)/2)) (e^(-(((π^2) τ)/4)))

So it is a simple matter to test various difference schemes by comparing with this known exact solution. It should be emphasized that this type of smooth initial data, which also joins continuously onto the boundary conditions, is rather atypical of option-pricing problems. Our purpose here is to simplify matters to get a general feel for the relative merits of explicit and implicit scheme.

Explicit Scheme

Here are the price- and time-steps used, together with an output that is the value of the parameter α=((Δτ)/(Δ (x^2))).

dx = 0.025; dtau = 0.00025; alpha = dtau/dx^2

The number of time-steps is 400, and there are 160 space-steps.

M=400; nminus = 80; nplus = 80;

Here we set initial and boundary conditions, and plot the latter:

initial=Table[N[Sin[Pi*(k-1-nminus)/nminus] ], {k,1, nminus+nplus+1}];
lower=Table[0, {m, 1, M+1}];
upper=Table[0, {m, 1, M+1}];

ListPlot[initial]

Here we define a function to solve the problem on the given grid, the output is a list of values of u for the last three time points (for now we shall just use the final value). Note that this function makes use of the Mathematica vector compilation methods. The initial conditions, and upper and lower boundary conditions, are supplied as vectors, initial, lower, upper, as rank-1 objects that are real - note the syntax {initial, _Real, 1} to denote this. For further details see section 2.6.15 of the Mathematica book. Note that you can specify a vector (rank 1), matrix (rank 2) and higher order objects. No statement about dimension is needed - just the rank of the "tensor".

ExplicitSolver =
	Compile[{{initial,_Real,1}, {lower,_Real,1}, {upper,_Real,1},alpha},
		Module[{wold=initial, wnew=initial, wvold=initial,m,k, tsize = Length[lower],xsize=Length[initial]},
			For[m=2, m <= tsize, m++,
				(wvold = wold; wold=wnew;
					For[k=2, k < xsize, k++,
					(wnew[[k]] = alpha (wold[[k-1]]+wold[[k+1]])+(1-2alpha) wold[[k]])];
						wnew[[1]] = lower[[m]];
						wnew[[xsize]] = upper[[m]])];
			{wvold,wold,wnew}]
		];

Now we apply this to get the solution at a time τ = 0.1 (the units are irrelevant for our analysis):

Timing[soln = ExplicitSolver[initial,lower,upper, alpha];]

Now we do the interpolation to supply a continuous function:

interpoldata = 
Table[{(k-nminus-1)*dx ,soln[[3,k]]}, {k, 1, nminus+nplus+1}];

ufunc = Interpolation[interpoldata, InterpolationOrder -> 3]

Now we plot the error in the answer:

Plot[ufunc[x]-Sin[((π x)/2)] (E^((1/4) (-(π^2)) 0.1`)),{x,-2,2},PlotPoints->50]

Finally we tabulate the numerical result, the exact result and the error in the numerical scheme. With 400 time-steps the error is manageably small.

samples = TableForm[Join[{{"x", "Explicit FD", "Exact", "Error"}},
Table[Map[PaddedForm[N[Chop[#1]],{5,6}]&,
N[{x, 
ufunc[x], 
Sin[Pi*x/2]*Exp[-(Pi^2 0.1)/4], 
ufunc[x]- Sin[Pi*x/2]*Exp[-(Pi^2 0.1)/4]},5]],
{x, -2, 2, 0.25}]]]

C++ Model of Explicit Diffusion Test Problem

Go and look at code explicitdiffusion2.cpp....compile and run it for this case. The following tells Mathematica where this code is on my demonstration laptop. It would need changing on any other system!

SetDirectory["C:\Documents and Settings\William Shaw\My Documents\LGS0809\cppexamples"]

FileNames[]

cppdata = ReadList["fdoutput.txt"];

ListPlot[cppdata]

Length[cppdata]

cppinterpoldata = 
Table[{(k-nminus-1)*dx ,cppdata[[k]]}, {k, 1, nminus+nplus+1}];

cppfunc = Interpolation[cppinterpoldata, InterpolationOrder -> 3]

Now we plot the error in the answer:

Plot[cppfunc[x]-Sin[((π x)/2)] (E^((1/4) (-(π^2)) 0.1`)),{x,-2,2},PlotPoints->50]

mmadata = Transpose[interpoldata][[2]];

mmacppdiff = mmadata-cppdata;

Max[Abs[mmacppdiff]]

These are identical to the precision we have given. They should be!

Unstable Explicit Scheme

Let's double the time step and take α to 0.8. Looks innocuous enough!

dx = 0.025; dtau = 0.0005; alpha = dtau/dx^2

The number of time-steps is 200, and there are 160 space-steps.

M=200; nminus = 80; nplus = 80;

M*dtau

Timing[soln = ExplicitSolver[initial,lower,upper, alpha];]

Now we do the interpolation to supply a continuous function:

interpoldata = 
Table[{(k-nminus-1)*dx ,soln[[3,k]]}, {k, 1, nminus+nplus+1}];

ufunc = Interpolation[interpoldata, InterpolationOrder -> 3]

Now we plot the error in the answer: (NOTE THE AXES!!)

Plot[ufunc[x]-Sin[((π x)/2)] (E^((1/4) (-(π^2)) 0.1`)),{x,-2,2},PlotPoints->50]

This is not a subtle matter. When you have instability, it is hard to miss, but you must do something to look for it. Now rerun the C++ code with M=200, N=80 and have a look at what happens.

cppdata = ReadList["fdoutput.txt", String];
cppdata[[Range[1,10]]]

Stable to Unstable Transition

Let's take a look at the growth of instability.

Contrast M=310,309,...300

In this demonstration I will make an animation to show the transition to instability. This will use the 
Manipulate function (new in Mathematica 6). If you have an older version then replace Manipulate by Do and you will get a sequence of graphics showing the same effect.

Manipulate[(M=310-k;
dx=0.025;dtau=(320/M)*0.00125/4;
nminus=80;nplus=80;alpha=dtau/dx^2;
soln = ExplicitSolver[initial,lower,upper, alpha];
interpoldata=Table[{(k-nminus-1)*dx,soln[[3,k]]},{k,1,nminus+nplus+1}];
ufunc=Interpolation[interpoldata,InterpolationOrder->3];
Plot[{ufunc[x],Sin[Pi*x/2]*Exp[-Pi^2 0.1/4]},{x,-2,2},PlotPoints->50, 
PlotRange -> {-2, 2}, PlotStyle -> {Thickness[0.001], Thickness[0.001]}, Epilog-> {Text[StringJoin["α = ",ToString[alpha]], {-1.5, 1.5}]}]), {k, 0, 7}]

Arrays in C++

This is an issue I am going to sidestep completely by using the public domain vector and matrix utilities created by the authors of Numerical Recipes in C++ [NRC++]. You might like to see the discussion of arrays in Section 1.2 of NRC++ under "Vectors and Matrices", as well as Appendix C of Joshi's book. Every C++ guru has probably written their own set. While I might have Mathematica guru status this does not remotely apply to C++ so I am going to work with the NRC++ set. Comments welcome on the relative merits of the schemes developed by NR, Joshi, Duffy and London etc.
[Now go to project with NR matrices and run that code - ExplicitDiffusion3]

SetDirectory["C:\Documents and Settings\William Shaw\My Documents\LGS0809\cppexamples"]

FileNames[]

cppdata2 = ReadList["fdoutput2.txt"];

mmacppdiff = mmadata-cppdata2;

Max[Abs[mmacppdiff]]

So this is fine too.