Abstract and Applied Analysis

Volume 2012, Article ID 120358, 20 pages

http://dx.doi.org/10.1155/2012/120358

## Double Discretization Difference Schemes for Partial Integrodifferential Option Pricing Jump Diffusion Models

Instituto de Matemática Multidisciplinar, Universitat Politècnica de València, Camino de Vera s/n, 46022 Valencia, Spain

Received 10 September 2012; Revised 7 November 2012; Accepted 7 November 2012

Academic Editor: Carlos Vazquez

Copyright © 2012 M.-C. Casabán et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

A new discretization strategy is introduced for the numerical solution of partial integrodifferential equations appearing in option pricing jump diffusion models. In order to consider the unknown behaviour of the solution in the unbounded part of the spatial domain, a double discretization is proposed. Stability, consistency, and positivity of the resulting explicit scheme are analyzed. Advantages of the method are illustrated with several examples.

#### 1. Introduction

Since empirical studies revealed that the normality of the log returns, as assumed by Black and Scholes, could not capture features like heavy tails and asymmetries observed in market-data log-returns densities [1], a number of models try to explain these empirical observations: stochastic volatility [2, 3], deterministic local volatility [4, 5], jump diffusion [6, 7], and infinite activity Lévy models [8–11]. The two last types of models, discussed in [12] and [13, chapters 14, 15] allow to calibrate the model to market price of options and reproduce a wide variety of implied volatility skews/smiles. These models are characterized by partial integrodifferential equations (PIDEs) that involve a second-order differential operator and a nonlocal integral term that requires specific treatment and presents additional difficulties.

In order to solve the PIDE problem numerically, Andersen and Andreasen [14] use an unconditionally stable ADI finite difference method and accelerate it using fast Fourier transform (FFT). In [15–17] wavelet methods are applied to infinite activity jump-diffusion models. Interesting analytic-numerical treatments for Lévy models have been introduced in [18–20]. The so-called COS method for pricing European options is presented in [18]. This is based on the knowledge of the characteristic function of the jump operator and the close relation of the characteristic function with the series coefficients of the Fourier-cosine expansion of the density function. In [19], an expansion of the characteristic function of local volatility models with Lévy jumps is developed. The authors in [20] derive an analytical formula for the price of European options for any model including local volatility and Poisson jump process by using Malliavin calculus techniques. Various authors apart of [14] used finite difference schemes for PIDEs in [21–27]. Discretization of the integral term leads to full matrices due to its nonlocal character. Dealing with the integral term several challenges arise, for instance, how to approximate the integral term and how to localizate a bounded computational domain, also the selection of the boundary conditions of the numerical domain and the problem of the double discretization of the differential and integral part of the PIDE.

Tavella and Randall in [26] used an implicit time discretization and propose a stationary fairly rapid convergent iterative method to solve the full matrix problem quoted above but without a careful numerical analysis. A generalization of this iterative method to price American options is proposed in [25].

In the outstanding paper [22] the authors propose an explicit-implicit finite difference scheme for solving parabolic PIDEs with possibly singular kernels when the random evolution of the underlying asset is driven by a time-inhomogeneous jump-diffusion process. The authors study stability and convergence of the proposed scheme as well as rates of convergence. However, they use backward or forward difference quotients of only first order depending on the sign of the coefficient of the convection term in order to avoid oscillations. An improvable issue of [22] is that in order to approximate the truncated integral term, they assume a particular behavior of the solution outside of the bounded numerical domain.

An efficient solution of PIDEs for the jump-diffusion Merton model is proposed in [24] with a very efficient treatment of the resulting dense linear system by using a circulant preconditioned conjugate gradient method. However, in [24], they only consider the particular case where the jump sizes have zero mean, . They also assume a particular behavior of the solution outside the bounded numerical domain.

Almendral and Oosterlee [28] present an implicit discretization of the PIDE jump-diffusion model on an uniform grid using finite differences, where a splitting technique combined with FFT is used to accelerate the dense matrix-vector product. The authors also assume a particular behaviour of the solution outside of the bounded numerical domain, in a similar way to [24].

In [21] a finite difference method for PIDE associated with the CGMY infinite activity Lévy model is treated. The equations are discretized in space by the collocation method and in time by an explicit backward differentiation formula. The integral part is transformed into a Volterra equation. After integration by parts and taking advantage of the vanishing derivative behaviour of the payoff function for large asset values, the authors are able to truncate properly the integral for the case of put and butterfly options.

In [27] the price of European and American options under PIDE Kou’s jump-diffusion model is solved using finite differences on nonuniform grids, and time stepping is performed using the implicit Rannacher scheme. The evaluation of the integral term is efficient from the computational cost point of view, assuming that the behaviour of the solution for large values of the underlying asset follows the asymptotic behaviour.

For the sake of clarity in the presentation we recall that in a jump-diffusion model, the modified stochastic differential equation (SDE) for the underlying asset is where is the underlying stock price, is the drift rate, is the volatility, is the increment of Gauss-Wiener process, and is the Poisson process. The random variable representing the jump amplitude is denoted by , and the expected relative jump size is denoted by . The jump intensity of the Poisson process is denoted by . Based on the SDE (1.1) the resulting PIDE for a contingent claim is given by [7, 14, 29]: where is the risk-free interest rate, the probability density of the jump amplitude is given by , and is the payoff function. Merton’s jump-diffusion model assumes that jump sizes are log-normally distributed with mean and standard deviation , that is,

In this paper we consider Merton’s jump-diffusion model for a vanilla call option with payoff function The aim of the paper is the construction and numerical analysis of an explicit finite difference numerical scheme of the PIDE (1.2)-(1.3), with a different treatment of the integral part from previously quoted authors. Instead of assuming long-term information about the solution we perform a full discretization of the integral part, involving the unknown function values in the numerical scheme, discriminating the finite truncation domain and the infinite remaining one. As a consequence, this strategy involves a double discretization with respect to the spatial variable. With respect to the time variable an explicit forward approximation is used. An account of the advantages of this explicit approach has been explained and applied in [30].

This paper is organized as follows. Section 2 deals with a transformation of variables in order to eliminate both the advection and the reaction terms of the PIDE (1.2). Then the integral part is split in two parts: a finite integral and an infinite one , and the last is again transformed into a finite integral. The separation point of the two split integrals is becoming a parameter that could be chosen according to the criteria used by [16, 22, 31]. A suitable choice of parameter is the one used by other authors when they truncate the numerical domain. For instance in [27], one takes ; in Section 6 we take . Section 3 deals with the construction of the numerical scheme and the selection of the numerical domain that always involves the difficulty of the consideration of the boundary conditions. For the case of a PIDE this issue is even more relevant because the values throughout all the unbounded integral domain are unknown. The spatial numerical domain is divided in two parts by the parameter : and . In the domain, the stepsize discretization is , consisting in equidistributed mesh points . The domain is transformed into the by transformation . In the transformed domain we consider a stepsize discretization and mesh points with . When the interval is reversed to the domain , the reversed mesh points become nonuniformly distributed. Hence, the numerical scheme for problem (2.9)-(2.10) is forward in time with time-step discretization . The approximation of is centered in of the unique parameter , and in the approximation of involves a nonuniform stepsize discretization depending on and the value of . The numerical approximation of the integrals is evaluated using trapezoidal quadrature rules with stepsizes and , respectively. The boundary conditions at the boundary of our numerical domain are as follows. At we assume that the solution is zero according to the vanilla call option problem. At our largest finite value considered we assume a linear behavior of the solution. This hypothesis has been previously used in [32].

In Section 4 sufficient conditions for stability and positivity of the numerical solution are given in terms of the three parameters and as well as of the parameter . Consistency of the scheme is treated in Section 5. Section 6 includes illustrative examples showing the possible advantages of our new discretization approach. Finally conclusions are shown in Section 7.

If is a vector in , we denote its infinite norm . Vector is said to be nonnegative if for all . Then we denote . For a matrix in , we denote by . Matrix is said to be nonnegative if for all , and we denote .

#### 2. Transformation of the Integrodifferential Problem

For the sake of convenience we introduce a transformation of variables to remove both the advection and the reaction terms of the PIDE problem (1.2)-(1.3).

Let us consider the transformation and note that problem (1.2)-(1.3) is transformed into the problem

In order to approximate the integral appearing in (2.2) and further discretization it is convenient the change of the variable

Let us denote

Taking , let us decompose Following [33, page 201] let us consider the substitution into , obtaining the expression Taking into account (2.4)–(2.8), the problem (2.2)-(2.3) can be written in the form

#### 3. Numerical Scheme Construction

In this section a difference scheme for problem (2.9)-(2.10) is constructed. With respect to the time variable, given with , let be the time-step discretization , and , , with integer. With respect to the spatial variable , given an arbitrary positive fixed , we construct a uniform grid in , with the spatial step discretization , with , being integer.

Note that the integral given by (2.6) requires the evaluation of the unknown at points . As this integral has been transformed into (2.8) over the interval ] for the variable , see (2.7), we consider a uniform mesh of into points, of the form , where is integer, , with . Taking into account (2.7), for the original variable in , one has and since , Thus the spatial domain is split into points from those only the first are equidistributed.

Let us denote For the approximation of we consider two types of finite differences: for the internal points of , and denoting , for the points lying in . Note that the discrete operator has different expressions depending on the ubication of , see (3.5) and (3.6).

From the previous approximations, for the internal points we have where and are approximations of the composite trapezoidal type of integrals appearing in (2.6), (2.8): Let us denote . The approximation takes the form where the first term for does not appear due to the null value of the limit of function given by (1.4) as tends to zero. On the other hand, considering the assumption that has asymptotic linear behaviour as () and using (1.4) it follows that the integrand of (3.9) verifies Consequently, the last term of related to vanishes, and one gets The numerical scheme (3.7)–(3.12) needs to incorporate the transformed initial condition and the boundary conditions for : and assuming linear behaviour of the solution for large values of the spatial variable at any time, we have and null integral term approximation . Hence, considering (3.7) for , one gets

For the sake of convenience in the study of stability we introduce a vector formulation of the scheme (3.7)–(3.15). Let us consider the vector in as and let be the tridiagonal matrix related to the differential part and defined by where

Let be the matrix in related to the integral part whose entries for each fixed in are defined by From the previous notation the scheme (3.7)–(3.15) can be written in the form

#### 4. Positive and Stability of the Numerical Solution

Dealing with prices of contracts modeled by PIDE, the solution must be nonnegative. In this section we show that numerical solution provided by scheme (3.7)–(3.15) is conditionally positive and stable.

We begin with the following result.

Lemma 4.1. *With previous notation, assume that stepsizes , in and , and in satisfy *(C1)*, *(C2)*. ** Then matrix given by (3.17) is nonnegative. *

*Proof. * From (3.18), for one has and for . On the other hand, for , we have that
Thus, under condition (C1), condition (4.1) holds true. With respect to the nonuniform grid, note that for , from (3.18) one gets that , , and . From (3.18) we also have that
In order to guarantee the nonnegativeness of the remaining entries of matrix , let us introduce the function
for . With this notation, we have that appearing in (3.18) satisfies
Note that
Taking into account that , with , then for , both numerator and denominator of (4.5) are positive. Thus
Thus is strictly increasing for , and its minimum is achieved at with the value
and condition (4.4) holds true under the condition
From condition (C2), properties (4.2) and (4.8) are satisfied and matrix is nonnegative.

Note that as matrix defined by (3.19) is always nonnegative, from Lemma 4.1 and (3.20) starting from nonnegative initial vector , the following result is established.

Theorem 4.2. * With the hypotheses and notation of Lemma 4.1, the solution of the scheme (3.7)–(3.15) is nonnegative if the initial values . *

The next result will be used below to guarantee stability.

Lemma 4.3. *Matrices and defined by (3.17), (3.18), and (3.19) satisfy the following. *(1)*Under conditions (C1) and (C2) of Lemma 4.1, . *(2)*, where . *

*Proof. * By Lemma 4.1, under hypotheses (C1) and (C2), all the entries of matrix are nonnegative. Thus,
We also have

Hence and from the definition of , it follows that . This proves part 1. From (3.19), for a fixed , with one gets where are the trapezoidal approximation rules with and points, approximating the integrals respectively.

Let be the maximum of in , given by Note that the integrand is increasing in the interval and decreasing for . In order to upper bound (4.12) we consider two cases. Firstly, let us assume that and let us denote to be the first integer with , such that From the properties of the lower Riemann sums, it follows that Taking into account the values and located just before and after , together with (4.17), from (4.12) it follows that In an analogous way, for the second situation where , one gets again (4.18). Finally, if is also true that Now we will upper bound (4.13). Let us denote with . Let be the maximum of : In this case we could distinguish the three possible situations ; , and upper bounding the lower Riemann sums relative to (4.13) one gets From (4.11), (4.18), and (4.22) together with the fact that , one gets for each value of , Taking into account that , for , it follows that and . Hence, from (4.23) one gets independently of the value of the size of matrix .

For the sake of clarity and as there are many definitions of stability in the literature we recall our concept of stability in the next definition.

*Definition 4.4. * Let be a numerical solution of the PIDE (2.9), (2.10) computed from the scheme (3.7)–(3.15) with stepsizes in , in , and in . Let be the corresponding vector form, that is, of (3.20). We say that is strongly uniformly stable, if
where is independent of , , , and .

If the property (4.25) is satisfied for appropriate relationships between the stepsizes , , and , then one says that the strong uniform stability is conditional.

Theorem 4.5. *With the previous notation, the numerical solution of the scheme (3.7)–(3.15) is strongly uniformly stable if one satisfies the condition together with
*

*Proof. *Note that scheme (3.7)–(3.15) is equivalent to the vector form scheme (3.20). Under condition (4.26), by Lemma 4.3 one gets, after taking norms in (3.20),
Hence, from (4.27), Bernouilli’s inequality, and , one gets
Thus the conditional strong uniform stability is established.

#### 5. Consistency

We say that a numerical difference scheme is consistent with a PIDE, if the exact theoretical solution of the PIDE approximates well to the exact solution of the difference scheme as the stepsize discretization tends to zero, [34, 35].

Let us write the scheme (3.7)–(3.12) in the form , where and let us write the PIDE (2.9) in the form where where and are given by (2.6)–(2.8).

Let us denote to be the value of the theoretical solution of PIDE (5.2). Let such that . We denote by the following expression the local truncation error : In order to prove the consistency, we must show that Assuming that is twice continuously partially differentiable with respect to and using Taylor’s expansions about one gets where

Let us assume that admits four times continuous partial derivatives with respect to , and let us denote In accordance with [34, page 101] let us explain the local consistency error of , see (2.6), by where in (5.11) means the second derivative with respect to the variable , see [33].

In an analogous way, let us explain the local consistency error of the unbounded integral , see (2.8), by Summarizing, one gets

Thus which proves the consistency of the scheme with the PIDE.

#### 6. Numerical Results

In the following examples the code was run on Matlab. The first example illustrates that stability conditions of Theorem 4.5 cannot be removed.

*Example 6.1. * Consider the vanilla call option problem (1.2)–(1.5) under Merton jump diffusion model with parameters , , , , , , and . Taking , , and , Figure 1 shows that when the stability conditions (4.26) are satisfied results are good , while if the stability conditions are broken results are unreliables.

The next example shows the robustness of our numerical scheme under changes of the jump intensity of the model.

*Example 6.2. * Taking the same parameters of Example 6.1 apart from and the stepsize discretizations , , and , Figure 2 shows the variation of the solution with parameter , where corresponds to the Black-Scholes case.

In the next example, the error is the difference between the numerical solution computed by (3.7)–(3.15) and (2.1) and the exact solution given by Merton’s formula [7]. Example 6.3 shows that the error of the numerical solution with fixed decreases with the uniform stepsize about the strike , while the error close to the truncation separation point remains stationary when decreases. This fact agrees with facts illustrated in [28, pages 15-16].

*Example 6.3. * Consider the vanilla call option problem (1.2)–(1.5) under Merton jump diffusion model with parameters , , , , , , and . For , , and , the Figure 3 shows the variation of the absolute error of the solution under changes of the stepsize .

The next Example 6.4 shows that the errors in the right boundary of the numerical domain when one uses finite difference schemes, quoted by [28], can be reduced with our double spatial discretization by decreasing the stepsize .

*Example 6.4. * Taking the problem of Example 6.3 with fixed , Figure 4 shows the error reduction of the numerical solution about the right boundary of the numerical domain when parameter decreases, while the error about the strike remains stationary.

#### 7. Conclusions

This work introduces a new discretization strategy for solving partial integrodifferential equations which involves the discretization of the unknown in the unbounded part of the integral. This fact increases the accuracy of the numerical solution in the boundary of the numerical domain as it is shown in Example 6.4.

#### Acknowledgment

This paper was supported by the Spanish M.E.Y.C. Grant DPI2010-20891-C02-01.

#### References

- J. Y. Campbell, A. W. Lo, and A. C. MacKinlay,
*The Econometrics of Financial Markets*, Princeton University Press, Princeton, NJ, USA, 1997. - S. Heston, “A closed-form solution for options with stochastic volatility with applications to bond and currency options,”
*Review of Financial Studies*, vol. 6, no. 2, pp. 327–343, 1993. View at Google Scholar - J. Hull and A. White, “The pricing of options with stochastic volatilities,”
*The Journal of Finance*, vol. 42, no. 2, pp. 281–300, 1987. View at Google Scholar - J. C. Cox and S. A. Ross, “The valuation of options for alternative stochastic processes,”
*Journal of Financial Economics*, vol. 3, no. 1-2, pp. 145–166, 1976. View at Google Scholar - B. Dupire, “Pricing with a smile,”
*Risk Magazine*, vol. 1, pp. 18–20, 1994. View at Google Scholar - S. G. Kou, “A jump diusion model for option pricing,”
*Management Science*, vol. 48, no. 8, pp. 1086–1101, 2002. View at Google Scholar - R. C. Merton, “Option pricing when the underlying stocks are discontinuous,”
*Journal of Financial Economics*, vol. 3, pp. 1125–2144, 1976. View at Google Scholar - O. E. Barndorff-Nielsen, “Processes of normal inverse Gaussian type,”
*Finance and Stochastics*, vol. 2, no. 1, pp. 41–68, 1998. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - E. Eberlein, “Application of generalized hyperbolic Lévy motions to finance,” in
*Lévy Processes|Theory and Applications*, O. Barndor-Nielsen, T. Mikosch, and S. Resnick, Eds., pp. 319–336, Birkhäuser, Boston, Mass, USA, 2001. View at Google Scholar · View at Zentralblatt MATH - I. Koponen, “Analytic approach to the problem of convergence of truncated Lévy ights towards the Gaussian stochastic process,”
*Physical Review E*, vol. 52, no. 1, pp. 1197–1199, 1995. View at Google Scholar - D. Madan and F. Milne, “Option pricing with variance gamma martingale components,”
*Mathematical Finance*, vol. 1, no. 4, pp. 39–55, 1991. View at Google Scholar - R. Cont and P. Tankov,
*Financial Modelling with Jump Processes*, Chapman & Hall/CRC Financial Mathematics Series, Chapman & Hall/CRC, Boca Raton, Fla, USA, 2004. - A. Pascucci,
*PDE and Martingale Methods in Option Pricing*, vol. 2 of*Bocconi & Springer Series*, Springer, Milan, Italy, 2011. View at Publisher · View at Google Scholar - L. Andersen and J. Andreasen, “Jump-diffusion processes: volatility smile fitting and numerical methods for option pricing,”
*Review of Derivatives Research*, vol. 4, no. 3, pp. 231–262, 2000. View at Google Scholar - A.-M. Matache, P.-A. Nitsche, and C. Schwab, “Wavelet Galerkin pricing of American options on Lévy driven assets,”
*Quantitative Finance*, vol. 5, no. 4, pp. 403–424, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - A.-M. Matache, T. von Petersdorff, and C. Schwab, “Fast deterministic pricing of options on Lévy driven assets,”
*Mathematical Modelling and Numerical Analysis*, vol. 38, no. 1, pp. 37–71, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - A.-M. Matache, C. Schwab, and T. P. Wihler, “Fast numerical solution of parabolic integrodifferential equations with applications in finance,”
*SIAM Journal on Scientific Computing*, vol. 27, no. 2, pp. 369–393, 2005. View at Publisher · View at Google Scholar - F. Fang and C. W. Oosterlee, “A novel pricing method for European options based on Fourier-cosine series expansions,”
*SIAM Journal on Scientific Computing*, vol. 31, no. 2, pp. 826–848, 2008/09. View at Publisher · View at Google Scholar - S. Pagliarani, A. Pascucci, and C. Riga, “Adjoint expansions in local Lévy models,” SSRN eLibrary, 2011.
- E. Benhamou, E. Gobet, and M. Miri, “Smart expansion and fast calibration for jump diffusions,”
*Finance and Stochastics*, vol. 13, no. 4, pp. 563–589, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - A. Almendral and C. W. Oosterlee, “Accurate evaluation of European and American options under the CGMY process,”
*SIAM Journal on Scientific Computing*, vol. 29, no. 1, pp. 93–117, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - R. Cont and E. Voltchkova, “A finite difference scheme for option pricing in jump diffusion and exponential Lévy models,”
*SIAM Journal on Numerical Analysis*, vol. 43, no. 4, pp. 1596–1626, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - Y. d'Halluin, P. A. Forsyth, and G. Labahn, “A penalty method for American options with jump diffusion processes,”
*Numerische Mathematik*, vol. 97, no. 2, pp. 321–352, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - E. W. Sachs and A. K. Strauss, “Efficient solution of a partial integro-differential equation in finance,”
*Applied Numerical Mathematics*, vol. 58, no. 11, pp. 1687–1703, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - S. Salmi and J. Toivanen, “An iterative method for pricing American options under jump-diffusion models,”
*Applied Numerical Mathematics*, vol. 61, no. 7, pp. 821–831, 2011. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - D. Tavella and C. Randall,
*Pricing Financial Instruments*, Wiley, New York, NY, USA, 2000. - J. Toivanen, “Numerical valuation of European and American options under Kou's jump-diffusion model,”
*SIAM Journal on Scientific Computing*, vol. 30, no. 4, pp. 1949–1970, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - A. Almendral and C. W. Oosterlee, “Numerical valuation of options with jumps in the underlying,”
*Applied Numerical Mathematics*, vol. 53, no. 1, pp. 1–18, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - M. Briani, C. La Chioma, and R. Natalini, “Convergence of numerical schemes for viscosity solutions to integro-differential degenerate parabolic problems arising in financial theory,”
*Numerische Mathematik*, vol. 98, no. 4, pp. 607–646, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - R. Company, L. Jódar, E. Ponsoda, and C. Ballester, “Numerical analysis and simulation of option pricing problems modeling illiquid markets,”
*Computers & Mathematics with Applications*, vol. 59, no. 8, pp. 2964–2975, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - R. Kangro and R. Nicolaides, “Far field boundary conditions for Black-Scholes equations,”
*SIAM Journal on Numerical Analysis*, vol. 38, no. 4, pp. 1357–1368, 2000. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - H. Windcli, P. A. Forsyth, and K. R. Vetzal, “Analysis of the stability of the stability of the linear boundary condition for the Black-Scholes equation,”
*Journal of Computational Finance*, vol. 8, no. 1, pp. 65–92, 2004. View at Google Scholar - P. J. Davis and P. Rabinowitz,
*Methods of Numerical Integration*, Computer Science and Applied Mathematics, Academic Press, New York, NY, USA, 2nd edition, 1984. - P. Linz,
*Analytical and Numerical Methods for Volterra Equations*, vol. 7, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, Pa, USA, 1985. View at Publisher · View at Google Scholar - G. D. Smith,
*Numerical Solution of Partial Differential Equations*, Clarendon Press, Oxford, UK, 3rd edition, 1985.