#### Daniel Duffy

##### C++ author, trainer

**Part I - Introduction and Initial Examples**

In this article we give an overview of the Boost C++ libraries (www.boost.org) and we discuss their applicability to certain problems in computational finance. There are approximately 90 libraries in Boost at the moment of writing and each one is focused on a particular aspect of software development. Part of the mission statement on www.boost.org reads:

*Boost provides free peer-reviewed portable C++ source libraries. We emphasize libraries that work well with the C++ Standard Library. Boost libraries are intended to be widely useful, and usable across a broad spectrum of applications. The Boost license encourages both commercial and non-commercial use.*

*We aim to establish "existing practice" and provide reference implementations so that Boost libraries are suitable for eventual standardization. Ten Boost libraries are already included in the C++ Standards Committee's Library Technical Report (TR1) and will be in the new C++0x Standard now being finalized. C++0x will also include several more Boost libraries in addition to those from TR1. More Boost libraries are proposed for TR2.*

The Boost libraries make use of templates and in this sense the generic programming model (GP) plays a more important role than the more traditional object-oriented programming (OOP) model. I have conveniently grouped the individual libraries into the following categories and I give some libraries in each category that I find useful:

*Higher-Order Functions*: define function types as objects. We can view these libraries as a generalization of and improvement on the use of function pointers in C. Libraries are: Function, Bind, Phoenix, Signals*Data Types*: data and number types: Tuple, Variant, Any, Rational, Integer, Interval. These are data types that improve the reliability or extend the functionality of existing C/C++ data types*Data Structures*: Complex structures that model geometrical and mathematical entities such as n-dimensional tensors, graphs and matrices (and matrix algebra)*Maths*: Univariate statistical distributions (about 25), special functions (e.g. gamma functions, Legendre polynomials), random number generation (e.g. Mersenne Twister), Accumulators (framework for statistical accumulators)*Text and string processing*: libraries that process, analyse and parse strings: StringAlgo, Tokenizer, Regex, Xpressive*Utilities, miscellaneous and other libraries*: many other useful libraries. We mention Smart Pointer, Thread, Date-Time, Serialization and Asio (portable networking)

In subsequent articles we shall discuss some specific Boost libraries and their applications in computational finance:

- Solving n-factor Black Scholes PDE with multi_array and finite differences
- Numerical linear algebra routines with uBLAS
- Why boost Function is so powerful
- Statistical distributions and special functions

The example we give is C++ code to calculate the value of pi (approximately 3.1415 for most quants). Here is the full and well-documented source code:

Code:

```
#include // Convenience header file
#include
#include // std::time
using namespace std;
#include
int main()
{
// A. Define the lagged Fibonacci and Uniform objects
boost::lagged_fibonacci607 rng;
rng.seed(static_cast (std::time(0)));
boost::uniform_real<> uni(0.0,1.0);
// B. Produce Uniform (0, 1)
boost::variate_generator >
uniRng(rng, uni);
// C. Choose the desired accuracy
cout << "How many darts to throw? "; long N; cin >> N;
// D. Thow the dart; does it fall in the circle or outside
// Start throwing darts and count where they land
long hits = 0;
double x, y, distance;
for (long n = 1; n <= N; ++n)
{
x = uniRng(); y = uniRng();
distance = sqrt(x*x + y*y);
if ( distance <=1)
{
hits++;
}
}
// E. Produce the answer
cout << "PI is: " << hits << ", " << 4.0 * double(hits) / double (N);
return 0;
}
```

**Part 2 - Functions, Functions, Functions**

**The Boost.Function Library**

- A Monte Carlo simulation that receives its random numbers from a variety of sources, for example Boost.Random (as discussed in Part I), the standard and unreliable
*rand()*function or a propriety in-house random-number generator, or an open-source library (for example, Quantlib) - An option pricing engine that interoperates with a range of numerical methods such as PDE/FDM, Monte Carlo and lattice models. The engine should have no knowledge of the internal workings of these numerical methods and it should use them solely by invoking their interfaces.

The steps when using Boost.Function in applications are:

- Specify signature the Boost.Function object (this entails its name, input arguments and return type).
- Instantiate the object from Step 1 by assigning it to a function object, class member function or C-style function
- Use the instantiated object in client code

*MyClient*that uses the exponential function

*exp(x)*, where x is a numeric type (typically a double, float or complex).

So, let’s get started. Since we are modelling all kinds of functions that accept a scalar value as input and return a scalar value we can specify this class of functions in Boost.Function as follows:

Code:

```
#include <boost/function.hpp>
typedef boost::function<double (double z)> ScalarFunction;
```

Code:

`typedef boost::function<double (double z1, double z2)> TwoParameterFunction;`

Code:

```
ScalarFunction myFunc;
TwoParameterFunction myFunc2;
```

Code:

```
// C function for exponential with 1 parameter
double Exp1d(double z)
{
return exp(z);
}
// C function for exponential with 2 parameters
double Exp2d(double z1, double z2)
{
return exp((z1+z2));
}
```

Code:

```
myFunc = Exp1d;
cout << myFunc(1.0) << endl;
myFunc2 = Exp2d;
cout << myFunc2(1.0, 0.5) << endl;
```

We now discuss a more complicated example that shows how powerful Boost.Function is. To this end, we create a class that has an embedded algorithm, part of which is customizable because it has a boost::function as private member. The interface is:

Code:

```
class MyClient
{ // A class with an embedded Boost.Function object. In
// this sense it is policy-free as it does not depend on
// a specific function form.
private:
// Function specification: output/input
ScalarFunction f;
public:
MyClient(const ScalarFunction& myFunc) { f = myFunc; }
double calculate(double z) const
{ // This simulates an algorithm
return 2.0 * f(z) + sin(z*z); // The actual function called here
}
// Change the implementation of the calculator
void SetFunction(const ScalarFunction& myFunc) { f = myFunc; }
};
```

Code:

```
MyClient myC(Exp1d);
cout << myC.calculate(x) << endl;
```

Code:

```
myC.SetFunction(AnotherFunc1d);
cout << myC1.calculate(x) << endl;
```

Code:

```
double AnotherFunc1d(double z)
{
return exp(z) * sin(z);
}
```

Code:

```
class ExponentialDecorator
{ // My own version of exponential with overload ()
// This is a function object
private:
double a;
public:
ExponentialDecorator(double param) { a = param;}
double operator () (double z) const
{ // This is the function call operator
return exp(a*z);
}
};
```

Code:

```
ExponentialDecorator myExp(20.0);
MyClient myC2(myExp);
cout << myC2.calculate(x) << endl;
ExponentialDecorator myExp2(1.0);
myC2.SetFunction(Exp1d);
cout << myC2.calculate(x) << endl;
```

The last example discusses the approximation of the exponential function by Padé rational functions which are quotients of two polynomials of degree q (the numerator) and p (the denominator). The corresponding template class has parameters representing the degrees q and p and is a function object:

Code:

```
template struct PadeApproximation
{
double operator ()(double z) const;
}
```

Code:

```
template <> struct PadeApproximation<1,0>
{ // Explicit Euler
double operator ()(double z) const
{
return 1.0 + z;
}
};
template <> struct PadeApproximation<0,1>
{ // Implicit Euler
double operator () (double z) const
{
return 1.0/(1.0 - z);
}
};
template <> struct PadeApproximation<1,1>
{ // Crank Nicolson
double operator () (double z) const
{
return (2.0 + z)/(2.0 - z);
}
};
```

Code:

```
double x = -0.09;
PadeApproximation<1,0> approxExplicitEuler;
PadeApproximation<1,0> approxImplicitEuler;
PadeApproximation<1,1> approxCrankNicolson;
MyClient myEE(approxExplicitEuler);
MyClient myIE(approxImplicitEuler);
MyClient myCN(approxCrankNicolson);
cout << myEE.calculate(x) << endl;
cout << myIE.calculate(x) << endl;
cout << myCN.calculate(x) << endl;
```

In the next installment we shall show how we have integrated Boost Function, multi_array and the ADE finite difference method to price two-factor and three-factor option PDE problems.

The full source code is on my web site http://www.datasimfinancial.com/forum/viewtopic.php?t=111

**Part 3 - Multi-dimensional Data Stuctures in Boost**

In this blog we introduce the Boost Multidimensional Array Library (MAL) that allows developers to define multidimensional arrays in C++ applications. The MultiArray concept defines an interface to hierarchically nested containers. In particular, we can define n-dimensional data structures using the template classes provided by the library and it provides operations to create multiarrays, access their elements, navigating in multiarrays as well as creating view of multiarray data. The underlying flexible memory model supports a number of layouts, thus making the library suitable for a wide range of applications. We postpone a discussion of applications until we have introduced MAL and given some examples. You can download the code for this blog from my site http://www.datasimfinancial.com/forum/viewtopic.php?t=111 and test the code for yourself.

We now give a first example of a three-dimensional hypercube structure. First, we assume that Boost has been installed on your system and that all necessary search paths have been defined (this aspect is system and compiler-dependent). Next, we need to include the following header file in our code:

Code:

`#include <boost/multi_array.hpp>`

Code:

```
// Define structure of tensor
const int dim = 3; // 3d matrix
typedef boost::multi_array<double, dim> Tensor;
Tensor tensor;
// Define the extents of each separate dimension
const int NT = 3;
const int NSIM = 4;
const int NDIM = 3;
tensor.resize(extents[NT][NSIM][NDIM]);
```

Code:

```
// Useful shorthand
Tensor tensor2(extents [NT][NSIM][NDIM]);
```

Code:

```
// Create a tensor using a Collections, in this case in Boost
array<Tensor::index, dim> extentsList = { { NT, NSIM, NDIM } };
Tensor tensor3 (extentsList);
```

Code:

```
// Define a 3d loop to initialise the data
typedef Tensor::index Index;
```

Code:

```
// Accessing elements using []
for (Index i = 0; i != NT; ++i)
{
for (Index j = 0; j != NSIM; ++j)
{
for (Index k = 0; k != NDIM; ++k)
{
tensor[i][j][k] = i + j + k;
}
}
}
```

Code:

```
// Accessing elements using ()
typedef boost::array Indices;
Indices indices;
for (indices[2]=0; indices[2] < NZ; indices[2]++)
{
for (indices[0]=0; indices[0] < NX; indices[0]++)
{
for (indices[1]=0; indices[1] < NY; indices[1]++)
{
tensor(indices) = 0.0;
}
}
}
```

As a last example, we wonder how a four-dimensional array can be created. In fact, it is not much more difficult than creating three-dimensional arrays. Here is an example:

Code:

```
// Typedef to make working with quads easier
const int quadDim = 4;
typedef boost::multi_array<double, quadDim> Quad;
// The dimensions
Quad::index dim1=3;
Quad::index dim2=4;
Quad::index dim3=2;
Quad::index dim4=2;
// Create a "quad" multi-array
Quad quad(boost::extents[dim1][dim2][dim3][dim4]);
// Fill the cube using the operator [index][index][index][index] syntax
for (Quad::index d1=0; d1<dim1; d1++)
{
for (Quad::index d2=0; d2<dim2; d2++)
{
for (Quad::index d3=0; d3<dim3; d3++)
{
for (Quad::index d4=0; d4<dim4; d4++)
{
quad[d1][d2][d3][d4]= 0.0;
}
}
}
}
```

**Part 4 - Unconditionally stable and second-order accurate explicit Finite Difference Method for Three Factor PDE**

**The ADE Scheme meets the Boost Libraries**

**Summary and Goals**

We now introduce the Alternating Direction Explicit ADE method to approximate the solution of the three-dimensional heat equation PDE in a bounded domain. ADE is a stable, second-order accurate finite difference scheme for general n-factor convection-diffusion equations (such as the Black Scholes PDE), it is easy to program and it is more efficient than conventional schemes such as Crank-Nicolson, ADI and Splitting.

The scope of this article is restricted to describing the ADE method, how we use it to solve a three-factor PDE and how we have implemented ADE using Boost. A future article will discuss the mathematical background to ADE and how we use it for approximating the solution of the n-factor Black-Scholes PDE.

The main goal of this note is to show how it is possible to solve high dimensional PDEs using a combination of the ADE method and Boost. In this sense we can model the problems using algebraic methods and we are closer to solving some problems associated with the

*curse of dimensionality*.

The focus in this blog is on using Boost libraries with ADE; we exclude a discussion of the numerical analysis and this topic will be discussed elsewhere.

**Background**

One of the numerical techniques used in derivatives pricing is the

*Finite Difference Method*(FDM) and it is used in Finance to approximate the solution of

*Partial Differential Equations*(PDE) that describe the time evolution of a derivative that is contingent on one or more underlying variables (Duffy 2006). The theory of FDM for one and two-factor problems is reasonably well developed but out interest in this note is to discuss the application of the

*Alternating Direction Explicit*(ADE) method to solving three factor derivatives PDE problems. In Duffy 2009 we discuss the applicability of ADE to one-factor option models where we motivate the method, discuss its mathematical and numerical foundations and we compare the scheme with other methods with regard to accuracy. We also give references to related work in Duffy 2009. Based on the results and successes with one-factor problems we have applied the method to multi-factor problems and in this blog we reduce the scope to the case of the three-dimensional heat equation as it contains many of the essential difficulties that we experience in multi-factor Black-Scholes PDEs.

Some of the advantages of using ADE for this class of PDE are:

- It is an unconditionally stable, second-order accurate and explicit finite difference scheme and this means that the method will give accurate results irrespective of the mesh sizes in the time and underlying (space) dimensions.
- ADE is efficient; in contrast to ADI and Splitting method we do not have to break a multi-factor PDE into a sequence of one-dimensional problems and in contrast to these latter methods we do not have to solve a tridiagonal matrix system at each time level and for each space dimension. The scheme is also very easy to design and to implement.
- We have applied ADE to multi-factor Black-Scholes PDEs using previously developed methods such as exponential fitting, the Janenko scheme for problems with mixed derivatives and the ability to transform the PDE in semi-infinite space to a PDE that is defined on a unit hypercube (Duffy 2009).
- It is easy to parallelise ADE by inserting compiler directives from the OpenMP library into the C++ code that implements the scheme, for example by using the parallel sections construct.
- ADE resolves the curse of dimensionality: this refers to the problems that we encounter when we try to design and implement finite difference schemes in three or more dimensions. Traditional finite difference schemes break down due to the complexity. The ADE scheme - in combination with Boost.MultAarray - resolves this curse to a large extent.

**The Model Problem**

In this section we discuss how to compute the solution of three-dimensional heat equation in a unit cube. The heat equation is one of the most important equations in mathematical physics and it has applications in many domains.

For this reason it is important to know how to find the solution (either in exact form or by numerical means). In this case we choose the finite difference method (FDM) to approximate the solution of the heat equation at equidistant mesh points in three-dimensional space. We assume that the reader knows what partial derivatives, partial differential equations (PDE) and divided differences are.

The PDE for the heat equation in three dimensions is given by:

(1)

[tex]\frac{\partial u}{\partial t} = \triangle u \equiv \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} + F(x,y,z,t), \quad x,y,z \in (0,1), \quad t > 0[/tex]There are an infinite number of solutions to (1). In order to specify a unique solution to (1) we must define some auxiliary conditions. First, we define the initial condition at t = 0:

(2)

[tex]u(x,y,z,0) = \varphi (x,y,z), \quad x,y,z \in (0,1)[/tex]In this case we state that the unknown solution u is known at t = 0. Second, we define zero

*Dirichlet boundary conditions*when the independent variables x, y and z are 0 and 1:

(3)

[tex]u(0,y,z)=u(1,y,z)=0, \quad y,z \in (0,1)[/tex][tex]u(x,0,z)=u(x,1,z)=0, \quad x,z \in (0,1)[/tex]

[tex]u(x,y,0)=u(x,y,1)=0, \quad x,y \in (0,1)[/tex]

We have now defined the initial boundary-value problem (1), (2), (3) defined on the unit cube and for all time t in the interval (0, T), where T > 0.

Turning to the finite difference solution of (1), (2), (3) we partition cube [0,1]^3 into three-dimensional mesh points.

To this end, we partition the interval [0,1] in each of the directions x, y and z into equal subintervals and we create three arrays as follows:

[tex]x_i=ih_1, \quad 0 \leq i \leq I[/tex]

[tex]y_j=jh_2, \quad 0 \leq j \leq J[/tex]

[tex]z_k=kh_3, \quad 0 \leq k \leq K[/tex]

where the mesh sizes are defined by:

[tex]h_1=1/I, \quad h_2=1/J, \quad h_3=1/K, I, J, K \geq 2[/tex] are integrers.

We also discretise the interval (0, T) into a number of equidistant

*mesh points*in time as follows:

[tex]\triangle t = T/NT, \quad t_n = n \triangle t, \quad 0 \leq n \leq NT[/tex]

Finally, we shall be working with mesh functions that are defined at the discrete mesh points in space and time:

[tex]v^n_{i,j,k} \sim v(x_i,y_j,z_k,t_n)[/tex]

In all further discussions we assume that all functions are mesh functions of this form.

**The ADE Scheme**

We are now ready to introduce the Alternating Direction Explicit (ADE) method (Duffy 2009). The scheme is stable, has second-order accuracy and it is explicit which means that we compute a solution without having to solve a matrix system at each time level. In general, we compute the final solution in a series of steps. First, we compute one component from the lower apex (0,0,0) to the upper apex (1,1,1) of the cube as follows. We call it the

*upward sweep:*

(4)

[tex]\frac{v^{n+1}_{i,j,k} - v^n_{i,j,k}}{\triangle t} = \frac{1}{h^2_1}(V^{n+1}_{i+1,j,k}-V^n_{i,j,k}-V^{n+1}_{i,j,k}+V^n_{i-1,j,k})[/tex][tex]+\frac{1}{h^2_1}(V^{n+1}_{i,j+1,k}-V^n_{i,j,k}-V^{n+1}_{i,j,k}+V^n_{i,j-1,k})[/tex]

[tex]+\frac{1}{h^2_3}(V^{n+1}_{i,j,k+1}-V^n_{i,j,k}-V^{n+1}_{i,j,k}+V^n_{i,j,k-1})[/tex]

[tex]1 \leq i \leq I-1, \quad 1 \leq j \leq J-1, \quad 1 \leq k \leq K-1, \quad n \geq 0[/tex]

The second component computes a solution from the upper apex to the lower apex. We call it the

*downward sweep*:

(5)

[tex]\frac{v^{n+1}_{i,j,k} - v^n_{i,j,k}}{\triangle t} = \frac{1}{h^2_1}(V^{n+1}_{i+1,j,k}-V^n_{i,j,k}-V^{n+1}_{i,j,k}+V^n_{i-1,j,k})[/tex][tex]+\frac{1}{h^2_1}(V^{n+1}_{i,j+1,k}-V^n_{i,j,k}-V^{n+1}_{i,j,k}+V^n_{i,j-1,k})[/tex]

[tex]+\frac{1}{h^2_3}(V^{n+1}_{i,j,k+1}-V^n_{i,j,k}-V^{n+1}_{i,j,k}+V^n_{i,j,k-1})[/tex]

[tex]1 \leq i \leq I-1, \quad 1 \leq j \leq J-1, \quad 1 \leq k \leq K-1, \quad n \geq 0[/tex] where [tex]h_{1}, h_{2}, h_{3}[/tex] are the mesh sizes in the x, y and z directions

Having computed these solutions (and this can be done in parallel) we arrive at the final solution as the average of the two components:

(6)

[tex]u^n_{i,j,k} = \frac{1}{2}(U^n_{i,j,k}+V^n_{i,j,k}) \quad 0 \leq i \leq I, \quad 0 \leq j \leq J, \quad 0 \leq k \leq K, \quad n \ge 0[/tex]In order to unambiguously define the solution (6) we define the discrete boundary conditions:

(7)

[tex]U^0_{i,j,k}= V^0_{i,j,k}=\varphi_{i,j,k} \equiv \varphi(x_i,y_j,z_k), \quad 1 \leq i \leq I-1,\quad 1 \leq j \leq J-1,\quad 1 \leq k \leq K-1[/tex]and the discrete boundary conditions for the component [tex]U^n_{i,j,k}[/tex] :

(8)

[tex]\left.\begin{aligned} U_{0,j,k}^n = U_{I,j,k}^n = 0, \; 0 \leq j \leq J, \; 0 \leq k \leq K\\ U_{i,0,k}^n = U_{i,J,k}^n = 0, \; 0 \leq i \leq I, \; 0 \leq k \leq K\\ U_{i,j,0}^n = U_{i,j,K}^n = 0, \; 0 \leq i \leq I, \; 0 \leq j \leq J \end{aligned}\right\rbrace n \geq 0[/tex]with the same boundary conditions holding for [tex]V^n_{i,j,k}[/tex]. Summarising, the scheme defined by equations (4) to (8) is an unambiguous specification of the finite difference method for the heat equation. We shall map this scheme using multi_array. To this end, we use some algebra to cast equations (4) and (5) to computable form:

(9)

[tex]U^{n+1}_{i,j,k}(1+\sigma)=U^{n}_{i,j,k}(1-\sigma)+\lambda_1(U^{n+1}_{i+1,j,k}+U^{n}_{i-1,j,k})+\lambda_2(U^{n+1}_{i,j+1,k}+U^{n}_{i,j-1,k})+\lambda_3(U^{n+1}_{i,j,k+1}+U^{n}_{i,j,k+1})+\Delta t F^{n+1}_{i,j,k}[/tex]and

(10)

[tex]V^{n+1}_{i,j,k}(1+\sigma)=V^{n}_{i,j,k}(1-\sigma)+\lambda_1(V^{n+1}_{i+1,j,k}+V^{n}_{i-1,j,k}) +\lambda_2(V^{n+1}_{i,j+1,k}+V^{n}_{i,j-1,k})+\lambda_3(V^{n+1}_{i,j,k+1}+V^{n}_{i,j,k+1})+\Delta t F^{n+1}_{i,j,k}[/tex]where

(11)

[tex]\lambda_j = \triangle t/h^2_j, \quad j=1,2,3,\quad \sigma= \sum_{j=1}^3 \lambda_j[/tex]**Design and Implementation in C++**

We now discuss how to set up the software system for solving the model PDE using the ADE scheme. The code that we show presently is based on a larger ongoing software project for the three-factor Black-Scholes equation. We have modified it for pedagogical reasons to make it readable.

The design approach is based on the strategy of decomposing the problem at hand into loosely coupled and autonomous (

*cohesive*) subsystems and classes. The main classes that we model are for:

- The coefficients of PDE (class
*ThreeFactorPde*). - The PDE domain in conjunction with its boundary and initial conditions (class
*ThreeFactorPdeDomain*). - The class that implements the ADE scheme (
*ThreeFactorADESolver*).

Code:

```
template class ThreeFactorHeatPde
{ // The parameter 'T' is usually instantiated to double in client code
public:
// Coefficients of the elliptic part of PDE
boost::function a11;
boost::function a22;
boost::function a33;
boost::function F; // RHS
ThreeFactorHeatPde() {}
};
```

Code:

```
template struct ThreeFactorPdeDomain
{
// 1. Domain
Range rx;
Range ry;
Range rz;
Range rt;
// 2. Boundary conditions, anticlockwise
// Specific boundaries in x
boost::function XLowerBC;
boost::function XUpperBC;
// Specific boundaries in y
boost::function YLowerBC;
boost::function YUpperBC;
// Specific boundaries in z
boost::function ZLowerBC;
boost::function ZUpperBC;
// 3. Initial condition
boost::function IC;
};
```

Code:

```
namespace ThreeFactorHeatImpl
{
// Domain information
double T;
double xMax, yMax, zMax; // Truncated domain
double diffusionX(double x, double y, double z, double t)
{
return 1.0;
}
double diffusionY(double x, double y, double z, double t)
{
return 1.0;
}
double diffusionZ(double x, double y, double z, double t)
{
return 1.0;
}
double H(double v)
{
return v*(1.0 - v);
}
double phi(double x, double y, double z)
{
return exp(-x*x) * exp(-y*y) * exp(-z*z);
}
double termInhomogeneous(double x, double y, double z, double t)
{
return (+6.0 - 4.0*(x*x+y*y+z*z)) * phi(x,y,z);
}
// Factories
void createThreeFactorHeatPde(ThreeFactorHeatPde& pde)
{
// Initialise all functions
pde.a11 = &ThreeFactorHeatImpl::diffusionX;
pde.a22 = &ThreeFactorHeatImpl::diffusionY;
pde.a33 = &ThreeFactorHeatImpl::diffusionZ;
pde.F = &ThreeFactorHeatImpl::termInhomogeneous;
}
double IC(double x, double y, double z)
{
return phi(x,y,z);
}
double BCXLower(double y, double z, double t)
{
return exp(-y*y) * exp(-z*z);
}
double BCXUpper(double y, double z, double t)
{
return exp(-xMax*xMax) * exp(-y*y) * exp(-z*z);
}
double BCYLower(double x, double z, double t)
{
return exp(-x*x) * exp(-z*z);
}
double BCYUpper(double x, double z, double t)
{
return exp(-x*x) * exp(-yMax*yMax) * exp(-z*z);
}
double BCZLower(double x, double y, double t)
{
return exp(-x*x) * exp(-y*y);
}
double BCZUpper(double x, double y, double t)
{
return exp(-x*x) * exp(-y*y) * exp(-zMax*zMax);
}
// Factory
void createThreeFactorPdeDomain(
ThreeFactorPdeDomain& pdeDomain)
{
pdeDomain.rx = Range(0.0, xMax);
pdeDomain.ry = Range(0.0, yMax);
pdeDomain.rz = Range(0.0, zMax);
pdeDomain.rt = Range(0.0, T);
pdeDomain.XLowerBC = &ThreeFactorHeatImpl::BCXLower;
pdeDomain.XUpperBC = &ThreeFactorHeatImpl::BCXUpper;
pdeDomain.YLowerBC = &ThreeFactorHeatImpl::BCYLower;
pdeDomain.YUpperBC = &ThreeFactorHeatImpl::BCYUpper;
pdeDomain.ZLowerBC = &ThreeFactorHeatImpl::BCZLower;
pdeDomain.ZUpperBC = &ThreeFactorHeatImpl::BCZUpper;
pdeDomain.IC = &ThreeFactorHeatImpl::IC;
}
} // End of namespace
```

Code:

```
class ThreeFactorADESolver
{ // Using ADE method for 3d heat equation PDE
private:
ThreeFactorHeatPde pde; // The pde coefficients
ThreeFactorPdeDomain pdeDomain; // Domain, BC, IC
// Data structures
BoostTensor* U; // upper sweep, n+1
BoostTensor* UOld; // upper sweep, n
BoostTensor* V; // lower sweep, n+1
BoostTensor* VOld; // lower sweep, n
public:
BoostTensor* MatNew; // averaged solution, level n+1
private:
// Mesh-related data
double hx, hy, hz, delta_k, // Step lengths in each direction
hx2, hy2, hz2, // 1/h^2
Sum, A11, A22, A33, // delta_k * (hx2 + hy2 + hz2)
G; // ~F // Right side of FDM for u_t = Lu + F
// Other variables
double tprev, tnow, T; // Time levels
long NX, NY,NZ,NT; // Number of subdivisions
double t1, t2, t3; // mesh values in each (x,y,z) direction
double coeffI, coeffII;
public:
// Mesh arrays in each direction
Vector xmesh;
Vector ymesh;
Vector zmesh;
Vector tmesh;
};
```

Code:

```
int const Dim = 3;
typedef boost::multi_array<double, Dim> BoostTensor;
typedef boost::array<BoostTensor::index, 3> BoostTensorIndices;
```

Code:

```
for (BoostTensor::index i = 1; i <= NX-1; i++)
{
for (BoostTensor::index j = 1; j <= NY-1; j++)
{
for (BoostTensor::index k = 1; k <= NZ-1; k++)
{
t1 = xmesh[i]; t2 = ymesh[j]; t3 = zmesh[k];
G = pde.F(t1,t2,t3,tnow) * delta_k;
(*U)[i][j][k] = (
(*UOld)[i][j][k]*coeffI
+ A11*((*UOld)[i+1][j][k] + (*U)[i-1][j][k])
+ A22*((*UOld)[i][j+1][k] + (*U)[i][j-1][k])
+ A33*((*UOld)[i][j][k+1] + (*U)[i][j][k-1])
+ G
)/coeffII;
}
}
}
for (BoostTensor::index i = NX-1; i >= 1; i--)
{
for (BoostTensor::index j = NY-1; j >= 1; j--)
{
for (BoostTensor::index k = NZ-1; k >= 1; k--)
{
t1 = xmesh[i]; t2 = ymesh[j]; t3 = zmesh[k];
G = pde.F(t1,t2,t3,tnow) * delta_k;
(*V)[i][j][k] = (
(*VOld)[i][j][k]*coeffI
+ A11*((*V)[i+1][j][k] + (*VOld)[i-1][j][k])
+ A22*((*V)[i][j+1][k] + (*VOld)[i][j-1][k])
+ A33*((*V)[i][j][k+1] + (*VOld)[i][j][k-1])
+ G
)/coeffII;
}
}
}
for (BoostTensor::index i = 0; i <= NX; i++)
{
for (BoostTensor::index j = 0; j <= NY; j++)
{
for (BoostTensor::index k = 0; k <= NZ; k++)
{
(*MatNew)[i][j][k] = 0.5*((*U)[i][j][k] + (*V)[i][j][k]); // == u
(*UOld)[i][j][k] = (*U)[i][j][k];
(*VOld)[i][j][k] = (*V)[i][j][k];
}
}
}
```

Code:

```
int main()
{
using namespace ThreeFactorHeatImpl;
xMax =1.0;
yMax = 1.0;
zMax = 1.0;
T = 2.0;
// Create pde and domain
ThreeFactorHeatPde pde;
createThreeFactorHeatPde(pde);
ThreeFactorPdeDomain pdeDomain;
createThreeFactorPdeDomain(pdeDomain);
cout << "Now creating solvern";
// FD scheme
long V = 100;
long NXX = V; long NYY = V; long NZZ = V; long NTT = V;
DatasimClock myClock;
myClock.start();
ThreeFactorADESolver solver(pde, pdeDomain, NXX, NYY, NZZ, NTT);
cout << "Now starting computationn";
solver.result();
myClock.stop();
cout << endl << "Running time " << myClock.duration() << endl;
BoostTensor exactTensor = CreateDiscreteFunction(&ExactSolution,
solver.xmesh, solver.ymesh, solver.zmesh,
T);
HotSpot maxError = MaxNorm(exactTensor, *(solver.MatNew));
// Get the spot where the largest error occurs
cout << endl << "Max error " << maxError.get<0>() << ", "
<< maxError.get<1>() << ", " << maxError.get<2>() << ", " << maxError.get<3>() << endl;
return 0;
}
```

Code:

```
double ExactSolution(double x, double y, double z, double t)
{
return exp(-x*x) * exp(-y*y) * exp(-z*z);
}
```

Code:

```
typedef tuple HotSpot;
HotSpot MaxNorm(const BoostTensor& tensor, const BoostTensor& tensor2)
{ // Compute largest element of m1 - m2 and state where it occurs
// For readability, define start and end indices for each dimension
BoostTensor::size_type size0 = tensor.shape()[0];
BoostTensor::index start0 = tensor.index_bases()[0];
BoostTensor::index end0 = start0 + size0;
BoostTensor::size_type size1 = tensor.shape()[1];
BoostTensor::index start1 = tensor.index_bases()[1];
BoostTensor::index end1 = start1 + size1;
BoostTensor::size_type size2 = tensor.shape()[2];
BoostTensor::index start2 = tensor.index_bases()[2];
BoostTensor::index end2 = start2 + size2;
double val = fabs(tensor[start0][start1][start2]
- tensor2[start0][start1][start2]);
double tmp;
HotSpot result;
for (BoostTensor::index row=start0; row val)
{
val = tmp;
// Update the hotspot
result.get<0>() = val;
result.get<1>() = row;
result.get<2>() = column;
result.get<3>() = layer;
}
}
}
}
return result;
}
```

**Summary**

This note builds on my previous two previous Quantnet.com blogs on the Boost libraries for higher-order functions and multi-dimensional arrays. We were particularly interested in showing how the ADE (Alternating Direction Explicit) method is applied to approximating the solution of the three-dimensional heat equation and implementing this scheme using Boost libraries.

The results from this feasibility study suggest our continuation of research into linear and nonlinear multi-factor models for equity and interest rate derivatives problem. We hope to report on these issues at a later stage.

**References**

- Duffy, D.J. 2006 Finite Difference methods in financial engineering John Wiley and Sons Chichester UK.
- Duffy, D.J. 2009 Unconditionally stable and second-order accurate explicit Finite Difference Schemes using Domain Transformation Part I. One-Factor Equity Problems. SSRN http://papers.ssrn.com/sol3/papers.cfm?abstract_id=1552926

**About the Author:**

Daniel Duffyis an author and trainer. His company Datasim specializes in methods and techniques for solving problems in quantitative finance. He is the author of Monte Carlo Frameworks: Building Customisable High-performance C++ Applications and Introduction to C++ for Financial Engineers: An Object-Oriented Approach. For more information on the author, see Quant Network's interview with Daniel Duffy

Daniel Duffy