Selected publications by Christina Christara

Supervised students' theses

Title: A high-order deferred correction method for the solution of free boundary problems using penalty iteration, with an application to American option pricing
Authors:Dawei Wang, Kirill Serkh, and Christina Christara
Date: May 2022, updated November 2022, March 2023
Pages: 41
Note: Submitted
Keywords: Green's functions; deferred correction; free boundary problems; immersed interface methods; finite differences; American options; Linear Complementarity Problem; penalty iteration; derivative jumps
File: HOFB.pdf

Abstract:

This paper presents a high-order deferred correction algorithm combined with penalty iteration for solving free and moving boundary problems, using a fourth-order finite difference method. Typically, when free boundary problems are solved on a fixed computational grid, the order of the solution is low due to the discontinuity in the solution at the free boundary, even if a high order method is used. Using a detailed error analysis, we observe that the order of convergence of the solution can be increased to fourth-order by solving successively corrected finite difference systems, where the corrections are derived from the previously computed lower order solutions. The corrections are applied solely to the right-hand side, and leave the finite difference matrix unchanged. The penalty iterations converge quickly given a good initial guess. We demonstrate the accuracy and efficiency of our algorithm using several examples. Numerical results show that our algorithm gives up to fourth-order convergence for both the solution and the free boundary location. We also test our algorithm on the challenging American put option pricing problem. Our algorithm gives the expected high-order convergence.

Title: Bilateral XVA pricing under stochastic default intensity: PDE modelling and computation
Authors:Yuwei Chen and Christina Christara
Date: June 2022
Pages: 32
Note:A later version is to be published in Applied Numerical Mathematics, vol 185, pages 236-259, March 2023. link to journal page
Keywords: Partial Differential Equations, Black-Scholes, Crank-Nicolson finite differences, mean reversion CIR process, stochastic default intensity, asymptotic approximation, nonlinear iteration, penalty iteration method.
File: sto_inten_xva_v3.pdf

Abstract: Adjusting derivative prices to take into account default risk has attracted the attention of several researchers and practitioners, especially after the 2007-2008 financial crisis. We derive a novel partial differential equation (PDE) model for derivative pricing including the adjustment for default risk, assuming that the default risk of one of the counterparties (the buyer) follows a Cox-Ingersoll-Ross (CIR) process, while the other party has constant default risk. The time-dependent PDE derived is of Black-Scholes type and involves two ``space'' variables, namely the asset price and the buyer default intensity, as well as a nonlinear source term. We formulate boundary conditions appropriate for the default intensity variable. The numerical solution of the PDE is based on standard finite differences, and a penalty-like iteration for handling the nonlinearity. We also develop and analyze a novel asymptotic approximation formula for the adjusted price of derivatives, resulting in a very efficient, accurate, and convenient for practitioners formula. We present numerical results that indicate stable second order convergence for the 2D PDE solution in terms of the discretization size. We compare the effectiveness of the 2D PDE and asymptotic approximations. We study the effect of various numerical and market parameters to the values of the adjusted prices and to the accuracy of the computed solutions.

Title: Penalty and Penalty-Like Methods for Nonlinear HJB PDEs
Authors:Christina Christara and Ruining Wu
Date: June 2022
Pages: 27
Note:A later version is published in Applied Mathematics and Computation, vol 425, pages 1-19, 2022. link to journal page
Keywords: Partial Differential Equations, Black-Scholes, nonlinear iteration, finite differences, Crank-Nicolson, control problem, Hamilton-Jacobi-Bellman (HJB) equation, transaction costs, stock borrowing fees, penalty methods, policy iteration.
File: hjb-penalty.pdf

Abstract: There are numerous financial problems that can be posed as optimal control problems, leading to Hamilton-Jacobi-Bellman or Hamilton-Jacobi-Bellman-Issacs equations. We reformulate these problems as nonlinear PDEs, involving max and/or min terms of the unknown function, and/or its first and second spatial derivatives. We suggest efficient numerical methods for handling the nonlinearity in the PDE through an adaptation of the discrete penalty method [Forsyth, Vetzal 2002] that gives rise to tridiagonal penalty matrices. We formulate a penalty-like method for the use with European exercise rights, and extend this to American exercise rights resulting in a double-penalty method. We also use our findings to improve the policy iteration algorithms described in [Forsyth, Labahn 2007]. Numerical results are provided showing clear second-order convergence, and where applicable, we prove the convergence of our algorithms.

Title: Penalty Methods for Bilateral XVA Pricing in European and American Contingent Claims by a PDE Model
Authors: Yuwei Chen and Christina C. Christara
Date: Nov 2019, Rev. Feb 2020
Pages: 22
Keywords: Partial Differential Equations, Black-Scholes, credit risk, default risk, credit valuation adjustment, call option, put option, forward contract, nonlinear iteration, finite differences, Crank-Nicolson, control problem, Hamilton-Jacobi-Bellman (HJB) equation.
Note: A version of this paper will appear in the
Journal of Computational Finance
File: ccc/xva_full.pdf

Abstract:
Taking into account default risk in the valuation of financial derivatives has become increasingly important, especially after the 2007-2008 financial crisis. Under some assumptions, the valuation of financial derivatives, including a value adjustment to account for default risk (the so-called XVA), gives rise to a nonlinear PDE \cite{Christoph Burgard and Mats Kjaer, 2011}. We propose numerical methods for handling the nonlinearity in this PDE, the most efficient of which are discrete penalty iteration methods. We first formulate a penalty iteration method for the case of European contingent claims, and study its convergence. We then extend the method to the case of American contingent claims, resulting in a double-penalty iteration. We also propose boundary conditions and their discretization for the XVA PDE problem in the cases of call and put options, as well as the case of a forward contract. Numerical results demonstrate the effectiveness of our methods.

Title: Analysis of quantization error in financial pricing via finite difference methods
Authors: Christina C. Christara and Nat Chun-Ho Leung
Date: July 2017, rev. March 2018
Pages: 27
Keywords: non-smooth initial conditions, option pricing, numerical solution, partial differential equation, convection-diffusion equations, Fourier analysis, finite difference methods
Note: A version of this paper is published in SIAM J. Numerical Analysis, 56:3, 2018, pp. 1731-1757
File: ccc/paper-error.pdf

Abstract:
In this paper, we study the error of a second order finite difference scheme for the one-dimensional convection-diffusion equation. We consider non-smooth initial conditions commonly encountered in financial pricing applications. For these initial conditions, we establish the explicit expression of the quantization error, which is loosely defined as the error of the numerical solution due to the placement of the point of non-smoothness on the numerical grid. Based on our analysis, we study the issue of optimal placement of such non-smoothness points on the grid, and the effect of smoothing operators on quantization errors.

Title: Partial Differential Equation pricing of contingent claims under stochastic correlation
Authors: Nat Chun-Ho Leung, Christina C. Christara and Duy-Minh Dang
Date: October 2017
Pages: 30
Keywords: stochastic correlation, option pricing, numerical solution, asymptotic solution, partial differential equation
Note: A version of this paper appeared in SIAM J. Scientific Computing, Vol 40, Issue 1, pages B1-B31 (2018).
File: ccc/paper_sto_corrL.pdf

Abstract:
In this paper, we study a partial differential equation (PDE) framework for option pricing where the underlying factors exhibit stochastic correlation, with an emphasis on computation. We derive a multi-dimensional time-dependent PDE for the corresponding pricing problem, and present a numerical PDE solution. We prove a stability result, and study numerical issues regarding the boundary conditions used. Moreover, we develop and analyze an asymptotic analytical approximation to the solution, leading to a novel computational asymptotic approach based on quadrature with a perturbed transition density. Numerical results are presented to verify second order convergence of the numerical PDE solution and to demonstrate its agreement The effect of certain problem parameters to the PDE solution, as well as to the asymptotic approximation solution, is also studied.

Title: Spread option pricing using ADI methods
Authors: Vida Heidarpour-Dehkordi and Christina C. Christara
Date: July 2017
Pages: 17
Keywords: Modified Craig-Sneyd, Alternating Direction Implicit method, two-dimensional Black-Scholes, American option, spread option, exchange option, analytical approximation, numerical PDE solution, penalty iteration
Note: Published in International J. Numerical Analysis and Modelling, 15:3, 2018. pp. 353-369.
File: ccc/spreadADI.pdf

Abstract:
Spread option contracts are becoming increasingly important, as they frequently arise in the energy derivative markets, e.g. exchange electricity for oil. In this paper, we study the pricing of European and American spread options. We consider the two-dimensional Black-Scholes Partial Differential Equation (PDE), use finite difference discretization in space and consider Crank-Nicolson (CN) and Modified Craig-Sneyd (MCS) Alternating Direction Implicit (ADI) methods for timestepping. In order to handle the early exercise feature arising in American options, we employ the discrete penalty iteration method, introduced and studied in Forsyth and Vetzal (2002), for one-dimensional PDEs discretized in time by the CN method. The main novelty of our work is the incorporation of the ADI method into the discrete penalty iteration method, in a highly efficient way, so that it can be used for two or higher-dimensional problems. The results from spread option pricing are compared with those obtained from the closed-form approximation formulae of Kirk (1995), Venkatramanan and Alexander (2011), Monte Carlo simulations, and the Brennan-Schwartz ADI Douglas-Rachford method, as implemented in MATLAB. In all spread option test cases we considered, including American ones, our ADI-MCS method, implemented on appropriate non-uniform grids, gives more accurate prices and Greeks than the MATLAB ADI method.

Title: PDE option pricing with variable correlations
Authors: Christina C. Christara and Nat Chun-Ho Leung
Date: May 2017
Pages: 17
Keywords: stochastic correlation, option pricing, regime switching, numerical solution, partial differential equation, system of PDEs, exchange option, basket option, Greeks.
Note: Proceedings of the International Conference on Computational and Informational Sciences and Engineering, honoring Professor Elias N. Houstis, University of Thessaly publications, September 2018, pp. 43-57.
File: ccc/paperRS.pdf

Abstract:
Modelling correlation between financial quantities is important in the accurate pricing of financial derivatives. In this paper, we introduce some stochasticity in correlation, by considering a regime-switching correlation model, in which the transition rates between regimes are given. We present a derivation of the associated Partial Differential Equation (PDE) problem. The problem involves a system of $\lsc$ PDEs, where $\lsc$ is the number of regimes. We formulate a finite difference method for the solution of the PDE system, and numerically demonstrate that it converges with second order. We study the effect of certain model parameters on the computed prices. We compare prices from this model, associated PDE and method with those from a stochastic correlation model, associated PDE and method in \cite{Emmerich06, Leung2017, LeungChristara2017} and discuss advantages and disadvantages.

Title: Option pricing in jump diffusion models with quadratic spline collocation
Authors: Christina C. Christara and Nat Chun-Ho Leung
Date: May 2014; revised May 2015
Pages: 15
Keywords: quadratic spline, collocation, finite difference, European option, American option, partial integro-differential equation, Merton model, Kou model, calculation of Greeks
Note: Published in Applied Mathematics and Computation, Vol 279, 10 April 2016, pp 28-42.
File: ccc/paper-jump.pdf,

Abstract:
In this paper, we develop a robust numerical method in pricing options, when the underlying asset follows a jump diffusion model. We demonstrate that, with the quadratic spline collocation method, the integral approximation in the pricing PIDE is intuitively simple, and comes down to the evaluation of the probabilistic moments of the jump density. The calculation of Greeks is also simple. When combined with a Picard iteration scheme, the pricing problem can be solved efficiently. We present the method and the numerical results from pricing European and American options with Merton's and Kou's models.

Title: An efficient numerical PDE approach for pricing foreign exchange interest rate hybrid derivatives
Authors: Duy Minh Dang, Christina C. Christara, Kenneth R. Jackson and Asif Lakhany
Date: March 2012; revised May 2013
Pages: 54
Keywords: Power-Reverse Dual-Currency (PRDC) swaps, Target Redemption (TARN), knockout, Partial Differential Equation (PDE), finite differences, non-uniform grids.
Note: Published in Journal of Computational Finance, Vol 18, Issue 4, 2015, pp 39-93.
File: prdcTARN.pdf, Link to the paper on the SSRN

Abstract:
We discuss efficient pricing methods via a Partial Differential Equation (PDE) approach for long-dated foreign exchange (FX) interest rate hybrids under a three-factor multi-currency pricing model with FX volatility skew. The emphasis of the paper is on Power-Reverse Dual-Currency (PRDC) swaps with popular exotic features, namely knockout and FX Target Redemption (FX-TARN). Challenges in pricing these derivatives via a PDE approach arise from the high-dimensionality of the model PDE, as well as from the complexities in handling the exotic features, especially in the case of the FX-TARN provision, due to its path-dependency. Our proposed PDE pricing framework for FX-TARN PRDC swaps is based on partitioning the pricing problem into several independent pricing sub-problems over each time period of the swap's tenor structure, with possible communication at the end of the time period. Each of these pricing sub-problems can be viewed as equivalent to a knockout PRDC swap with a known time-dependent barrier, and requires a solution of the model PDE, which, in our case, is a time-dependent parabolic PDE in three space dimensions. Finite difference schemes on non-uniform grids are used for the spatial discretization of the model PDE, and the Alternating Direction Implicit (ADI) timestepping methods are employed for its time discretization. Numerical examples illustrating the convergence properties and efficiency of the numerical methods are provided.

Title: Graphics processing unit pricing of exotic cross-currency interest rate derivatives with a foreign exchange volatility skew model,
Authors:Duy Minh Dang, Christina C. Christara and Kenneth R. Jackson
Date: February 2011
Pages: 17
Keywords: Power Reverse Dual Currency (PRDC) Swaps, Bermudan Cancelable, Partial Differential Equation (PDE), Alternating Direction Implicit (ADI), Finite Differences, Graphics Processing Units (GPUs), GPU Clusters, MPI, Parallel Computing
Note: Published in Journal of Concurrency and Computation: Practice and Experience (JCCPE) Volume 26, Issue 9, 2014, pp.~1609--1625.
File: prdc.gpu.bc.mod.pdf, Link to the paper on the SSRN

Abstract:

We present a Graphics Processing Unit (GPU) parallelization of the computation of the price of exotic cross-currency interest rate derivatives via a Partial Differential Equation (PDE) approach. In particular, we focus on the GPU-based parallel pricing of long-dated foreign exchange (FX) interest rate hybrids, namely Power Reverse Dual Currency (PRDC) swaps with Bermudan cancelable features. We consider a three-factor pricing model with FX volatility skew which results in a time-dependent parabolic PDE in three spatial dimensions. for the spatial discretization of the PDE, and the Alternating Direction Implicit (ADI) technique is employed for the time discretization. We then exploit the parallel architectural features of GPUs together with the Compute Unified Device Architecture (CUDA) framework to design and implement an efficient parallel algorithm for pricing PRDC swaps. Over each period of the tenor structure, we divide the pricing of a Bermudan cancelable PRDC swap into two independent pricing subproblems, each of which can efficiently be solved on a GPU via a parallelization of the ADI timestepping technique. Numerical results indicate that GPUs can provide significant increase in performance over CPUs when pricing PRDC swaps. An analysis of the impact of the FX skew on such derivatives is provided.

Title:A Highly Efficient Implementation on GPU Clusters of PDE-Based Pricing Methods for Path-Dependent Foreign Exchange Interest Rate Derivatives
Authors:Duy Minh Dang, Christina C. Christara and Kenneth R. Jackson
Date: Submitted March 2013, revised April 2013.
Pages: 17
Keywords: Power Reverse Dual Currency (PRDC) Swaps, Bermudan Cancelable, Partial Differential Equation (PDE), Alternating Direction Implicit (ADI), Finite Differences, Graphics Processing Units (GPUs), GPU Clusters, MPI, Parallel Computing
Note: Published in Springer Lecture Notes in Computer Science, Volume 7975, 2013, pp.~107--126, Proceedings of The 13th International Conference on Computational Science and Its Applications (ICCSA 2013), Ho Chi Minh City, Vietnam, 24--27 June 2013.
File: dang-ccc-krj-iccsa13.pdf, Link to the paper on the SSRN

Abstract:

We present a highly efficient parallelization of the computation of the price of exotic cross-currency interest rate derivatives with path-dependent features via a Partial Differential Equation (PDE) approach. In particular, we focus on the parallel pricing on Graphics Processing Unit (GPU) clusters of long-dated foreign exchange (FX) interest rate derivatives, namely Power-Reverse Dual-Currency (PRDC) swaps with FX Target Redemption (FX-TARN) features under a three-factor model. Challenges in pricing these derivatives via a PDE approach arise from the high-dimensionality of the model PDE, as well as from path-dependency of the FX-TARN feature. The PDE pricing framework for FX-TARN PRDC swaps is based on partitioning the pricing problem into several independent pricing sub-problems over each time period of the swap's tenor structure, with possible communication at the end of the time period. Finite difference methods on non-uniform grids are used for the spatial discretization of the PDE, and the Alternating Direction Implicit (ADI) technique is employed for the time discretization. Our implementation of the pricing procedure on a GPU cluster involves (i) efficiently solving each independent sub-problems on a GPU via a parallelization of the ADI timestepping technique, and (ii) utilizing MPI for the communication between pricing processes at the end of the time period of the swap's tenor structure. Numerical results showing the efficiency of the parallel methods are provided.

Title:An Efficient GPU-Based Parallel Algorithm for Pricing Multi-Asset American Options
Authors:Duy Minh Dang, Christina C. Christara and Kenneth R. Jackson
Date:March 2011
Pages: 14
Keywords: American Option, Multi-Asset, Penalty Method, Alternating Direction Implicit Approximate Factorization (ADI-AF), Graphics Processing Units, GPUs, Parallel Computing, Finite Difference
Note: A revised version of this paper (with a slightly different title) appeared in the Journal of Concurrency and Computation: Practice and Experience (JCCPE) special issue on high-performance computational finance, Vol. 24, No. 8, June 2012, pp. 849-866.
File: Link to the paper on the SSRN

Abstract:

We develop highly-efficient parallel Partial Differential Equation (PDE) based pricing methods on Graphics Processing Units (GPUs) for multi-asset American options. Our pricing approach is built upon a combination of a discrete penalty approach for the linear complementarity problem arising due to the free boundary and a GPU-based parallel Alternating Direction Implicit Approximate Factorization technique with finite differences on uniform grids for the solution of the linear algebraic system arising from each penalty iteration. A timestep size selector implemented efficiently on GPUs is used to further increase the efficiency of the methods. We demonstrate the efficiency and accuracy of the parallel numerical methods by pricing American options written on three assets.

Title:Adaptive and high-order methods for valuing American options
Authors:Christina C. Christara and Duy Minh Dang
Date:December 2009
Pages: 24
Keywords: adaptive mesh selection, error equidistribution, quadratic splines, collocation, finite differences, European option, American option, penalty method
Note: A later version of this paper is published in the Journal of Computational Finance, Vol 14, Issue 4, 2011, pages 73-113.
File: ccc/americanOption.pdf

Abstract:
We develop space-time adaptive and high-order methods for valuing American options using a partial differential equation approach. The linear complementarity problem arising due to the free boundary is handled by a penalty method. Both finite difference and finite element methods are considered for the space discretization of the PDE, while classical finite differences, such as Crank-Nicolson, are used for the time discretization. The high-order discretization in space is based on an optimal finite element collocation method, the main computational requirements of which are the solution of one tridiagonal linear system at each time step, while the resulting errors at the grid points and midpoints of the space partition are fourth-order. To control the space error, we use adaptive gridpoint distribution based on an error equidistribution principle. A time stepsize selector is used to further increase the efficiency of the methods. Numerical examples show that our methods converge fast and provide highly accurate options prices, Greeks, and early exercise boundaries.

Title:A Parallel Implementation on GPUs of ADI Finite Difference Methods for Parabolic PDEs with Applications in Finance
Author:Duy Minh Dang, Christina C. Christara and Kenneth R. Jackson
Date:March 2010; revised December 2010; published 2011
Pages: 21
Keywords: Alternating Direction Implicit, ADI, Partial Differential Equation, PDE, Graphics Processing Units, GPUs, parallel computing, finite difference, multi-asset options, multi-currency swaps
Note: Appeared in the Canadian Applied Mathematics Quarterly (CAMQ), Vol. 17, No. 4, 2009, pp. 627-660, which was published in 2011, for a belated 30-year anniversary issue of CAMQ.
File: Link to the paper on the SSRN

Abstract:

We study a parallel implementation on a Graphics Processing Unit (GPU) of Alternating Direction Implicit (ADI) time-discretization methods for solving time-dependent parabolic Partial Differential Equations (PDEs) in three spatial dimensions with mixed spatial derivatives in a variety of applications in computational finance. Finite differences on uniform grids are used for the spatial discretization of the PDEs. As examples, we apply the GPU-based parallel methods to price European rainbow and European basket options, each written on three assets. Numerical results showing the efficiency of the parallel methods are provided.

Title:A PDE Pricing Framework for Cross-Currency Interest Rate Derivatives
Author:Duy Minh Dang, Christina C. Christara, Kenneth R. Jackson and Asif Lakhany
Date:November 2009
Pages: 10
Keywords: Power Reverse Dual Currency Swaps, Bermudan Cancelable, Partial Differential Equation, PDE, Alternating Direction Implicit, ADI, Generalized Minimal Residual, GMRES, Fast Fourier Transform, FFT
Note: Submitted to the "2010 Workshop on Computational Finance and Business Intelligence", International Conference In Computational Science (ICCS), Amsterdam, May 31 to June 2, 2010.
File: Link to the paper on the SSRN

Abstract:
We propose a general framework for efficient pricing via a Partial Differential Equation (PDE) approach of cross-currency interest rate derivatives under the Hull-White model. In particular, we focus on pricing long-dated foreign exchange (FX) interest rate hybrids, namely Power Reverse Dual Currency (PRDC) swaps with Bermudan cancelable features. We formulate the problem in terms of three correlated processes that incorporate FX skew via a local volatility function. This formulation results in a parabolic PDE in three spatial dimensions. Finite difference methods are used for the spatial discretization of the PDE. The Crank-Nicolson method and the Alternating Direction Implicit (ADI) method are considered for the time discretization. In the former case, the preconditioned Generalized Minimal Residual (GMRES) method is employed for the solution of the resulting block banded linear system at each time step, with the preconditioner solved by Fast Fourier Transform (FFT) techniques. Numerical examples illustrating the convergence properties of the methods are provided.

Title:Pricing Multi-Asset American Options on Graphics Processing Units Using a PDE Approach
Authors:Duy Minh Dang, Christina C. Christara and Kenneth R. Jackson
Date:September 2010
Pages: 8
Keywords: American Option, Multi-Asset, Penalty Method, Alternating Direction Implicit Approximate Factorization (ADI-AF), Graphics Processing Units, GPUs, Parallel Computing, Finite Difference
Note: Proceedings of the Third Workshop on High Performance Computational Finance at SC10, The International Conference for High Performance Computing, Networking, Storage, and Analysis, sponsored by ACM/IEEE, New Orleans, USA, 13-19 November 2010.
File: DOI:10.1109/WHPCF.2010.5671831
ccc/american.gpu.ieee.pdf, link to pdf file in ieeexplore

Abstract:

We develop highly efficient parallel pricing methods on Graphics Processing Units (GPUs) for multi-asset American options via a Partial Differential Equation (PDE) approach. The linear complementarity problem arising due to the free boundary is handled by a penalty method. Finite difference methods on uniform grids are considered for the space discretization of the PDE, while classical finite differences, such as Crank-Nicolson, are used for the time discretization. The discrete nonlinear penalized equations at each timestep are solved using a penalty iteration. A GPU-based parallel Alternating Direction Implicit Approximate Factorization technique is employed for the solution of the linear algebraic system arising from each penalty iteration. We demonstrate the efficiency and accuracy of the parallel numerical methods by pricing American options written on three assets.

Title: Quartic spline collocation for second-order boundary value problems
Authors:Christina C. Christara and Guohong Liu
Date:September 2009
Pages: 8
Keywords: sixth order convergence, quartic splines, second-order BVP, collocation
Note: Proceedings of the 9th HERCMA Conference on Computer Mathematics and Applications, Athens University of Economics and Business, September 23-26, 2009, Athens, Greece, pp. 1-8.
File: ccc/paper-qrtbvp2.pdf

Abstract:
Collocation methods based on quartic splines are presented for second-order two-point boundary value problems. In order to obtain a uniquely solvable linear system for the degrees of freedom of the quartic spline collocation approximation, in addition to the boundary conditions specified by the problem, extra boundary or near-boundary conditions are introduced. Non-optimal (fourth-order) and optimal (sixth-order) quartic-spline methods are considered. The theoretical behavior of the collocation methods is verified by numerical experiments. The extension of the methods to two-dimensional problems is briefly considered.

Title: Quadratic spline collocation for one-dimensional parabolic partial differential equations
Authors: Christina C. Christara, Tong Chen and Dang Duy Minh
Date: 2009
Pages: 43
Note: Published in Numerical Algorithms.
File:DOI doi:10.1007/s11075-009-9317-9
Journal page: http://www.springerlink.com/content/v1g8565702373446/
Earlier version: http://www.cs.toronto.edu/pub/reports/na/ccc/parab.ps.gz
Keywords: quadratic splines, collocation, parabolic PDEs, Crank-Nicolson, stability, optimal order of convergence, American options.

Abstract:

New methods for solving general linear parabolic partial differential equations (PDEs) in one space dimension are developed. The methods combine quadratic-spline collocation for the space discretization and classical finite differences, such as Crank-Nicolson, for the time discretization. The main computational requirements of the most efficient method are the solution of one tridiagonal linear system at each time step, while the resulting errors at the gridpoints and midpoints of the space partition are fourth order. The stability and convergence properties of some of the new methods are analyzed for a model problem. Numerical results demonstrate the stability and accuracy of the methods. Adaptive mesh techniques are introduced in the space dimension, and the resulting method is applied to the American put option pricing problem, giving very competitive results.

Title:Quartic spline collocation for fourth-order boundary value problems
Authors: Christina C. Christara, Ying Zhu and Jingrui Zhang
Date:September 2008
Pages: 6
Note: Proceedings of the 2008 Numerical Analysis conference, September 1-5, 2008, Kalamata, Greece, pp. 62-67.
File: ccc/paper-qrtbvp4.pdf

Abstract:
We consider the numerical solution of linear fourth-order boundary value problems (BVPs) in one and two dimensions, by methods based on quartic splines and the collocation methodology. The discretization error is sixth order on the gridpoints and midpoints of a uniform partition and fifth order globally in the uniform norm. For the linear systems arising from the discretization of certain biharmonic problems by quartic spline collocation, we develop fast solvers based on Fourier transforms, leading to asymptotically almost optimal solution techniques.

Title:A High-Performance Method for the Biharmonic Dirichlet Problem on Rectangles
Authors:Christina C. Christaramand Jingrui Zhang
Date: July 2007
Pages: 3
Note: Proceedings of the 2007 International Conference On Preconditioning Techniques for Large Sparse Matrix Problems in Scientific and Industrial Applications, July 9-12, 2007 Météopole, Toulouse, France, pp. 35-37.
File: Please contact ccc SYMBOL_AT cs DOT toronto DOT edu

Abstract:
We propose a fast solver for the linear system resulting from the application of a sixth-order Bi-Quartic Spline Collocation method to the biharmonic Dirichlet problem on a rectangle. The fast solver is based on Fast Fourier Transforms (FFTs) and preconditioned GMRES (PGMRES) and has complexity $O(n^2\log(n))$ on an $n \times n$ uniform partition. The FFTs are applied to two auxiliary problems with different boundary conditions on the two vertical sides of the domain, while the PGMRES is applied to a problem related to the two vertical boundaries. We show that the number of PGMRES iterations required to reduce the relative residual to a certain tolerance is independent of the gridsize $n$. Numerical experiments verify the effectiveness of the solver.

Title: Optimal Quadratic and Cubic Spline Collocation on Nonuniform Partitions
Authors: Christina C. Christara, and Kit Sun Ng
Date: 2006
Pages: 31
Note:Appeared in Computing, 76:3-4, 2006, pp. 227-257.
File:SCnonuni.pdf, DOI
Keywords: spline collocation; second-order two-point boundary value problem; error bounds; optimal order of convergence; spline interpolation.

Abstract:

We develop optimal quadratic and cubic spline collocation methods for solving linear second-order two-point boundary value problems on non-uniform partitions. To develop optimal non-uniform partition methods, we use a mapping function from uniform to non-uniform partitions and develop expansions of the error at the non-uniform collocation points of some appropriately defined spline interpolants. The existence and uniqueness of the spline collocation approximations are shown, under some conditions. Optimal global and local orders of convergence of the spline collocation approximations and derivatives are derived, similar to those of the respective methods for uniform partitions. Numerical results on a variety of problems, including a boundary layer problem, and a non-linear problem, verify the optimal convergence of the methods, even under more relaxed conditions than those assumed by theory.

Title: Adaptive Techniques for Spline Collocation
Authors: Christina C. Christara, and Kit Sun Ng
Date: 2006
Pages: 19
Note:Appeared in Computing, 76:3-4, 2006, pp. 259-277.
File:adapt.pdf, DOI
Keywords: spline collocation; second-order two-point boundary value problem; error bounds; optimal order of convergence; adaptive grid; grid size estimator; error estimator; spline interpolation.

Abstract:

We integrate optimal quadratic and cubic spline collocation methods for second-order two-point boundary value problems with adaptive grid techniques, and grid size and error estimators. Some adaptive grid techniques are based on the construction of a mapping function that maps uniform to non-uniform points, placed appropriately to minimize a certain norm of the error. One adaptive grid technique for cubic spline collocation is mapping-free and resembles the technique used in COLSYS (COLNEW). Numerical results on a variety of problems, including problems with boundary or interior layers, and singular perturbation problems indicate that, for most problems, the cubic spline collocation method requires less computational effort for the same error tolerance, and has equally reliable error estimators, when compared to Hermite piecewise cubic collocation. Comparison results with quadratic spline collocation are also presented.

Title: GIA-induced secular variations in the Earth's long wavelength gravity field: Influence of 3-D viscosity variations (short communication),
Authors: Konstantin Latychev, Jerry X. Mitrovica, Mark E. Tamisiea, Jeroen Tromp, Christina C. Christara and Robert Moucha
Date: November 2005
Pages: 6
Note:Published in Earth and Planetary Science Letters, Vol. 240, Issue 2, Pages 322-327
File:DOI link , journal page
Keywords: geopotential harmonics, glacial isostatic adjustment, 3-D structure

Abstract:

Predictions of present day secular variations in the Earth's long wavelength geopotential driven by glacial isostatic adjustment (GIA) have previously been analyzed to infer the radial profile of mantle viscosity and to constrain ongoing cryospheric mass balance. These predictions have been based on spherically symmetric Earth models. We explore the impact of lateral variations in mantle viscosity using a new finite-volume formulation for computing the response of 3-D Maxwell viscoelastic Earth models. The geometry of the viscosity field is constrained from seismic-to-mographic images of mantle structure, while the amplitude of the lateral viscosity variations is tuned by a free parameter in the modeling. We focus on the zonal J_l harmonics for degrees l = 2, ..., 8 and demonstrate that large-scale lateral viscosity variations of two to three orders of magnitude have a modest, 5-10%, impact on predictions of J_2. In contrast, predictions of higher degree harmonics show a much greater sensitivity to lateral variation in viscosity structure. We conclude that future analyses of secular trends (for degree l > 2) estimated from ongoing (GRACE, CHAMP) satellite missions must incorporate GIA predictions based on 3-D viscoelastic Earth models.

Title: Glacial isostatic adjustment on 3-D Earth models: a finite-volume formulation
Authors: Konstantin Latychev, Jerry X. Mitrovica, Jeroen Tromp, Mark E. Tamisiea, Dimitri Komatitsch, Christina C. Christara
Date: May 2005
Pages: 24
Note:Published in Geophysical Journal International, Vol. 161, Issue 2, Pages 421-444
File:DOI link , journal page
Keywords: crystal motions, finite volumes, glacial isostatic adjustment, 3-D viscoelastic Earth models, parallel computation, domain decomposition

Abstract:

We describe and present results from a finite-volume (FV) parallel computer code for forward modelling the Maxwell viscoelastic response of a 3-D, self-gravitating, elastically compressible Earth to an arbitrary surface load. We implement a conservative, control volume discretization of the governing equations using a tetrahedral grid in Cartesian geometry and a low-order, linear interpolation. The basic starting grid honours all major radial discontinuities in the Preliminary Reference Earth Model (PREM), and the models are permitted arbitrary spatial variations in viscosity and elastic parameters. These variations may be either continuous or discontinuous at a set of grid nodes forming a 3-D surface within the (regional or global) modelling domain. In the second part of the paper, we adopt the FV methodology and a spherically symmetric Earth model to generate a suite of predictions sampling a broad class of glacial isostatic adjustment (GIA) data types (3-D crustal motions, long-wavelength gravity anomalies). These calculations, based on either a simple disc load history or a global Late Pleistocene ice load reconstruction (ICE-3G), are benchmarked against predictions generated using the traditional normal-mode approach to GIA. The detailed comparison provides a guide for future analyses (e.g. what grid resolution is required to obtain a specific accuracy?) and it indicates that discrepancies in predictions of 3-D crustal velocities less than 0.1 mm yr^{-1} are generally obtainable for global grids with ~3 x 10^6 nodes; however, grids of higher resolution are required to predict large-amplitude (>1 cm yr^{-1}) radial velocities in zones of peak post-glacial uplift (e.g. James bay) to the same level of absolute accuracy. We conclude the paper with a first application of the new formulation to a 3-D problem. Specifically, we consider the impact of mantle viscosity heterogeneity on predictions of present-day 3-D crustal motions in North America. In these tests, the lateral viscosity variation is constructed, with suitable scaling, from tomographic images of seismic S-wave heterogeneity, and it is characterized by approximately 2 orders of magnitude (peak-to-peak) lateral variations within the lower mantle and 1 order of magnitude variations in the bulk of the upper mantle (below the asthenosphere). We find that the introduction of 3-D viscosity structure has a profound impact on horizontal velocities; indeed, the magnitude of the perturbation (of order 1 mm yr^{-1}) is as large as the prediction generated from the underlying (1-D) radial reference model and it far exceeds observational uncertainties currently being obtained from space-geodetic surveying. The relative impact of lateral viscosity variations on predicted radial motions is significantly smaller.

Title: Optimal Quadratic Spline Collocation on Non-Uniform Partitions
Authors: Christina C. Christara, and Kit Sun Ng
Date: June 2003, revised 2004
Pages: 28
Note:A revised version appeared in Computing.
File: ccc/QSCnonuni.ps.gz
Keywords: spline collocation; second-order two-point boundary value problem; error bounds; optimal order of convergence; adaptive grid; grid size estimator; error estimator; spline interpolation.

Abstract:

Quadratic and Cubic Spline Collocation (QSC and CSC) methods of optimal orders of convergence have been developed for the solution of linear second-order two-point Boundary Value Problems (BVPs) discretized on uniform partitions. In this paper, we extend the optimal QSC methods to non-uniform partitions, and in a companion paper we do the same for CSC. To develop optimal non-uniform partition QSC methods, we use a mapping function from uniform to non-uniform partitions and develop expansions of the error at the non-uniform collocation points of some appropriately defined spline interpolants. The existence and uniqueness of the QSC approximations are shown, under some conditions. The $j$th derivative of the QSC approximation is $O(h^{3-j})$ globally, and $O(h^{4-j})$ locally on certain points, for $j \ge 0$.

Title: Optimal Cubic Spline Collocation on Non-Uniform Partitions
Authors: Christina C. Christara, and Kit Sun Ng
Date: June 2003, revised 2004
Pages: 17
Note:A revised version appeared in Computing.
File: ccc/CSCnonuni.ps.gz
Keywords: spline collocation; second-order two-point boundary value problem; error bounds; optimal order of convergence; adaptive grid; grid size estimator; error estimator; spline interpolation.

Abstract:

In a recent paper, non-uniform partition Quadratic Spline Collocation (QSC) methods of optimal order of convergence were developed for second-order two-point Boundary Value Problems (BVPs). In this paper, we develop optimal non-uniform Cubic Spline Collocation (CSC) methods. The formulation and analysis of the methods are based, as in the QSC case, on a mapping function from uniform to non-uniform partitions. However, the CSC method is turned into a mapping-free method. The $j$th derivative of the CSC approximation is $O(h^{4-j})$ globally, for $j \ge 0$, and $O(h^{5-j})$ locally on certain points, for $j > 0$. The non-uniform partition CSC method is integrated with adaptive grid techniques, and grid size and error estimators. One adaptive grid technique is based on the construction of a mapping function. Another is mapping-free and resembels the technique used in COLSYS (COLNEW) \cite{AMR95}. Numerical results on a variety of problems, including problems with boundary or interior layers, and singular perturbation problems indicate that, for most problems, the CSC method require less computational effort for the same error tolerance, and has equally reliable error estimators, when compared to Hermite piecewise cubic collocation (HPCC).

Title: Quadratic Spline Galerkin Method for the Shallow Water Equations
Authors:Anita T. Layton, Christina C. Christara, and Kenneth R. Jackson
Date: November 2002
Pages: 30
Note: Mathematics and Computers in Simulation, 71:3, 2006, pp. 175-186
File: sweqsg02.ps.gz
Keywords: numerical weather prediction; finite element; semi-Lagrangian semi-implicit; Rossby stability; staggered grids
Abstract:
Currently in most global meteorological applications, the spectral transform method or low-order finite difference/finite element methods are used. The spectral transform method, which yields high-order approximations, requires Legendre transforms. The Legendre transforms have a computational complexity of $\mathcal O(N^3)$, where $N$ is the number of subintervals in one dimension, and thus render the spectral transform method unscalable. In this study, we present an alternative numerical method for solving the shallow water equations (SWEs) on a sphere in spherical coordinates. In this implementation, the SWEs are discretized in time using the two-level semi-Lagrangian semi-implicit method, and in space on staggered grids using the quadratic spline Galerkin method. We show that, when applied to a simplified version of the SWEs, the method yields a neutrally stable solution for the meteorologically significant Rossby waves. Moreover, we demonstrate that the Helmholtz equation arising from the discretization and solution of the SWEs should be derived algebraically rather than analytically, in order for the method to be stable with respect to the Rossby waves. These results are verified numerically using Boyd's equatorial wave equations \cite{boyd1} with initial conditions chosen to generate a soliton.

Title: Optimal Quadratic Spline Collocation Methods for the Shallow Water Equations
Authors:Anita T. Layton, Christina C. Christara, and Kenneth R. Jackson
Date: November 2002
Pages: 38
Note: Mathematics and Computers in Simulation, 71:3, 2006, pp. 187-205
File: sweqsc02.ps.gz
Keywords: numerical weather prediction; finite element; semi-Lagrangian semi-implicit; Rossby stability; staggered grids
Abstract:
In this study, we present numerical methods, based on the optimal quadratic spline collocation (OQSC) methods, for solving the shallow water equations (SWEs) in spherical coordinates. A quadratic spline collocation method approximates the solution of a differential problem by a quadratic spline. In the standard formulation, the quadratic spline is computed by making the residual of the differential equations zero at a set of collocation points; the resulting error is second order, while the error associated with quadratic spline interpolation is fourth order locally at certain points and third order globally. The OQSC methods generate approximations of the same order as quadratic spline interpolation. In the one-step OQSC method, the discrete differential operators are perturbed to eliminate low-order error terms, and a high-order approximation is computed using the perturbed operators. In the two-step OQSC method, a second-order approximation is generated first, using the standard formulation, and then a high-order approximation is computed in a second phase by perturbing the right sides of the equations appropriately.
In this implementation, the SWEs are discretized in time using the semi-Lagrangian semi-implicit scheme, which allows large timesteps while maintaining numerical stability, and in space using the OQSC methods. The resulting methods are efficient and yield stable and accurate representation of the meteorologically important Rossby waves. Moreover, by adopting the Arakawa C-type grid, the methods also faithfully capture the group velocity of inertia-gravity waves.

Title: Quadratic Spline Collocation Revisited: Extension to Systems of Elliptic PDEs,
Authors:Christina C. Christara and Kit Sun Ng
Date: July 2002
Pages: 32
Note: Depart. of Computer Science, University of Toronto, Technical Report 318/2001, rev. July 2002.
File: ccc/sys.ps.gz

Abstract:
We revisit the optimal Quadratic Spline Collocation (QSC) methods for single elliptic PDEs developed in \cite{Chri94} and prove some new results, including a discrete maximum principle. We extend QSC to systems of two coupled linear second-order PDEs in two dimensions. The system of PDEs is treated as one entity; no decoupling is applied.
We study the matrix properties of the linear system arising from the discretization of systems of two PDEs by QSC. We give sufficient conditions under which the QSC linear system is uniquely solvable and develop analytic formulae for the eigenvalues and eigenvectors of the QSC matrix from $2 \times 2$ Helmholtz operators with constant coefficients.
Numerical results demonstrate that the QSC methods are optimal, that is, fourth order locally at certain points and third order globally, even for some problems that do not satisfy the conditions of the analysis. The QSC methods are compared to conventional second-order discretization methods and are shown to produce smaller approximation errors in the same computation time, while they achieve the same accuracy in less time.
The formulation of the optimal two-step QSC method for a $2 \times 2$ system of PDEs can be extended in a straightforward way to a $n \times n$ system of PDEs.

Title: Fast Fourier Transform Solvers and Preconditioners for Quadratic Spline Collocation
Authors: Christina C. Christara and Kit Sun Ng
Date: November 2000
Pages: 33
Note: Depart. of Computer Science, University of Toronto, Technical Report 317/2000. A revised version to appear in BIT, 42:4 (December 2002).
File: ccc/precond.pdf, ccc/precond.ps.gz
DOI: 10.1023/A:1021944218806
Journal: PDF file (fulltext), page
Abstract:
Quadratic Spline Collocation (QSC) methods of optimal order of convergence have been recently developed for the solution of elliptic Partial Differential Equations (PDEs). In this paper, linear solvers based on Fast Fourier Transforms (FFT) are developed for the solution of the QSC equations. The complexity of the FFT solvers is O(N^2 log N), where N is the gridsize in one dimension. These direct solvers can handle PDEs with coefficients in one variable or constant, and Dirichlet, Neumann and periodic boundary conditions. General variable coefficient PDEs are handled by preconditioned iterative solvers. The preconditioner is the QSC matrix arising from a constant coefficient PDE. The convergence analysis of the preconditioner is presented. It is shown that, under certain conditions, the convergence rate is independent of the gridsize. The preconditioner is solved by FFT techniques, and integrated with one-step or acceleration methods, giving rise to asymptotically almost optimal linear solvers, with complexity O(N^2 log N). Numerical experiments verify the effectiveness of the solvers and preconditioners, even on problems more general than the analysis assumes. The development and analysis of FFT solvers and preconditioners is extended to QSC equations corresponding to systems of elliptic PDEs.

Title: An Efficient Transposition Algorithm for Distributed Memory Computers
Authors: C. Christara, X. Ding and K. R. Jackson
Date: October 1999
Pages: 19
Note: Proceedings of the 13th Annual International Symposium on High Performance Computing Systems and Applications (HPCS'99), Queen's University, Kingston, Ontario, Canada, 13-16 June 1999.
File: transpose.99b.ps.gz
Abstract:
Data transposition is required in many numerical applications. When implemented on a distributed-memory computer, data transposition requires all-to-all communication, a time consuming operation. The Direct Exchange algorithm, commonly used for this task, is inefficient if the number of processors is large. We investigate a series of more sophisticated techniques: the Ring Exchange, Mesh Exchange and Cube Exchange algorithms. These data transposition schemes were incorporated into a parallel solver for the shallow-water equations. We compare the performance of these schemes with that of the Direct Exchange Algorithm and the MPI all-to-all communication routine, MPI_AllToAll. The numerical experiments were performed on a Cray T3E computer with 512 processors and on an ethernet-connected cluster of 36 Sun workstations. Both the analysis and the numerical results indicate that the more sophisticated Mesh and Cube Exchange algorithms perform better than either the simpler well-known Direct and Ring Exchange schemes or the MPI_AllToAll routine. We also generalize the Mesh and Cube Exchange algorithms to a d-dimensional mesh algorithm, which can be viewed as a generalization of the standard hypercube data transposition algorithm.

Title: Multigrid and Multilevel Methods for Quadratic Spline Collocation
Authors: Christina C. Christara and Barry Smith
Date: February 1997
Pages: 21
Note: Appeared in BIT 37:4 (1997), pp. 781-803.
File: multign.pdf, multign.ps.gz
DOI: 10.1007/BF02510352
Journal: PDF file (fulltext), page
Abstract:
Multigrid methods are developed and analyzed for quadratic spline collocation equations arising from the discretization of one-dimensional second-order differential equations. The rate of convergence of the two-grid method integrated with a damped Richardson relaxation scheme as smoother is shown to be faster than 1/2, independently of the step-size. The additive multilevel versions of the algorithms are also analyzed. The development of quadratic spline collocation multigrid methods is extended to two-dimensional elliptic partial differential equations. Multigrid methods for quadratic spline collocation methods are not straightforward: because the basis functions used with quadratic spline collocation are not nodal basis functions, thus the design of efficient restriction and extension operators is nontrivial. Experimental results, with V-cycle and full multigrid, indicate that suitably chosen multigrid iteration is a very efficient solver for the quadratic spline collocation equations.

Title: Parallel solvers for spline collocation equations
Author: C. C. Christara
Date: January 1995
Pages: 35
Note: Appeared in Advances in Engineering Software, 27:1/2, 1996, pp. 71-89
File: ccc/civil.ps.gz
Abstract:
A variety of solvers for the spline collocation equations arising from the discretisation of elliptic Partial Differential Equations (PDEs) are considered. The convergence properties of semi-iterative and Krylov subspace acceleration methods applied to the system of spline collocation equations are studied. The preconditioners tested include incomplete factorisation and SSOR for both the natural and multicolour orderings, domain decomposition based on Schur complement methods with nonoverlapping subdomains, or Schwarz methods with overlapping subdomains, and multigrid methods.

The parallelisation of some of the above iterative methods is studied and their advantages and disadvantages discussed. The communication requirements of the methods are discussed when the methods are implemented on distributed memory machines. Results which show that spline collocation methods are very competitive for the solution of PDEs are presented.

Title: Quadratic spline collocation methods for elliptic partial differential equations
Author: C. C. Christara
Date: March 1994 (replaces Tech. Rep. DCS-135, 1990)
Pages: 34
Note: A slightly revised version of this paper appeared in BIT, 34:1, 1994, pp. 33-61
File: ccc/cs-90-135.ps.gz
DOI: 10.1007/BF01935015
Journal: PDF file (fulltext), page
Abstract:
We consider Quadratic Spline Collocation (QSC) methods for linear second order elliptic Partial Differential Equations (PDEs). The standard formulation of these methods leads to non-optimal approximations. In order to derive optimal QSC approximations, high order perturbations of the PDE problem are generated. These perturbations can be applied either to the PDE problem operators or to the right sides, thus leading to two different formulations of optimal QSC methods. The convergence properties of the QSC methods are studied. Optimal $O(h sup 3-j )$ global error estimates for the $j$-th partial derivative are obtained for a certain class of problems. Moreover, $O(h sup 4-j )$ error bounds for the $j$-th partial derivative are obtained at certain sets of points. Results from numerical experiments verify the theoretical behaviour of the QSC methods. Performance results also show that the QSC methods are very effective from the computational point of view. The QSC methods have been implemented efficiently on parallel machines.

Title: High performance computing of elliptic partial differential equations with spline collocation
Author: C. C. Christara
Date: June 1994
Pages: 13
Note: Appeared in Proceedings of the 1994 Computational Structures Technology conference (CST94), August 30 - September 1, 1994, Athens, Greece, sponsored by Civil-Comp, Vol. A, pp. 23-34
File: ccc/highperf.ps.gz
Abstract:
We consider a variety of solvers for the spline collocation equations arising from the discretisation of elliptic Partial Differential Equations. We study the convergence properties of semi-iterative and conjugate gradient acceleration methods applied to the system of spline collocation equations. The preconditioners tested include incomplete factorisation and SSOR for both the natural and multicolour orderings, domain decomposition based on Schur complement methods with nonoverlapping subdomains, or Schwarz methods with overlapping subdomains, and multigrid methods. We study the parallelisation of some of the above iterative methods and discuss their advantages and disadvantages. The communication requirements of the methods are discussed when the methods are implemented on distributed memory machines. Results which show that spline collocation methods are very competitive for the solution of BVPs are presented.

Title: Parallel computation of partial differential equations on distributed memory machines
Author: C. C. Christara
Date: June 1992
Pages: 13
Note: Appeared in Advances in Computer Methods for Partial Differential Equations - VII, pp. 142-149 R. Vichnevetsky, D. Knight, G. Richter eds.
File: ccc/paral.ps.gz
Abstract:
We review a number of parallel methods for the solution of elliptic Boundary Value Problems (BVPs). The methods considered are of coarse or medium grain parallelism and fixed grid discretisation. The cases where the degree of parallelism grows linearly with the size of the problem in all its dimensions or in fewer dimensions are considered. We also consider the parallelism in both the discretisation of the BVP phase and the solution of the resulting linear system phase. We discuss the assignment of data to processors and the communication requirements of the methods, when implemented on distributed memory machines. Based on a general model for a distributed memory machine, we develop a framework in order to classify the communication requirements of the methods considered. We study the communication complexity of the methods on a general distributed machine model, and present performance results on the iPSC/2 hypercube.

Title: Multicolour orderings and iterative methods for elliptic equations
Author: C. C. Christara
Date: June 1991
Pages: 11
Note: Appeared in Proceedings of the Supercomputing Symposium '91 (SS '91), June 3-5, 1991, Fredericton, NB, Canada, pp. 221-230
File: ccc/colour.ps.gz
Abstract:
In this paper we study the performance and parallel implementation of iterative methods applied to linear systems arising from the discretisation of linear elliptic Boundary Value Problems (BVPs) which are ordered according to multicolour orderings. We view the discretisation of a BVP as a process of generating a set of equations between stencils. These equations, when ordered according to some indexing of the stencil, result in a sparse linear system. We discuss the effect of indexing of the stencils on the sparsity pattern of the linear system, on the convergence of iterative methods applied to the system, as well as on the parallel implementation of these methods. We discuss which stencil indexings are preferable in order to minimise communication in message-passing parallel architectures.

Title: Domain decomposition and incomplete factorisation methods for partial differential equations
Author: C. C. Christara
Date: May 1991
Pages: 11
Note: Appeared in Proceedings of the sixth Distributed Memory Computing Conference (DMCC6), April 28 - May 1, 1991, Portland, OR, U.S.A., pp. 369-377
File: ccc/iccg.ps.gz
Abstract:
In this paper we develop and study a method which tries to combine the merits of Domain Decomposition (DD) and Incomplete Cholesky preconditioned Conjugate Gradient method (ICCG) for the parallel solution of linear elliptic Partial Differential Equations (PDEs) on rectangular domains. We first discretise the PDE problem, using Spline Collocation, a method of Finite Element type based on smooth splines. This gives rise to a sparse linear system of equations. The ICCG method provides us with a very efficient, but not straightforward parallelisable linear solver for such systems. On the other hand, DD methods are very effective for elliptic PDEs. A combination of DD and ICCG methods, in which the subdomain solves are carried out with ICCG, leads to efficient and highly parallelisable solvers. We implement this hybrid DD\-ICCG method on a hypercube, discuss its parallel efficiency, and show results from experiments on configurations with up to 32 processors. We apply a totally local communication scheme and discuss its performance on the iPSC/2 hypercube. A similar approach can be used with other PDE discretisation methods.

Title: Schur complement preconditioned conjugate gradient methods for spline collocation equations
Author: C. C. Christara
Date: May 1991
Pages: 14
Note: An earlier version of this paper appeared in Computer architecture news, 18:3, 1990, pp. 108-120, and Proceedings of the 1990 International Conference on Supercomputing, June 1990, Amsterdam, the Netherlands, sponsored by ACM
File: ccc/scmpcg.ps.gz
Abstract:
We are interested in the efficient solution of linear second order Partial Differential Equation (PDE) problems on rectangular domains. The PDE discretisation scheme used is of Finite Element type and is based on quadratic splines and the collocation methodology. We integrate the Quadratic Spline Collocation (QSC) discretisation scheme with a Domain Decomposition (DD) technique. We develop DD motivated orderings of the QSC equations and unknowns and apply the Preconditioned Conjugate Gradient (PCG) method for the solution of the Schur Complement (SC) system. Our experiments show that the SC-PCG-QSC method in its sequential mode is very efficient compared to standard direct band solvers for the QSC equations. We have implemented the SC-PCG-QSC method on the iPSC/2 hypercube and present performance evaluation results for up to 32 processors configurations. We discuss a type of nearest neighbour communication scheme, in which the amount of data transfer per processor does not grow with the number of processors. The estimated efficiencies are at the order of 90%.

Title: Conjugate gradient methods for spline collocation equations
Author: C. C. Christara
Date: May 1991
Pages: 9
Note: An earlier version of this paper appeared in Proceedings of the fifth Distributed Memory Computing Conference (DMCC5), April 8-12, 1990, Charleston, SC, U.S.A., pp. 550-558
File: ccc/pcgcoef.ps.gz
Abstract:
We study the parallel computation of linear second order elliptic Partial Differential Equation (PDE) problems in rectangular domains. We discuss the application of Conjugate Gradient (CG) and Preconditioned Conjugate Gradient (PCG) methods to the linear system arising from the discretisation of such problems using quadratic splines and the collocation discretisation methodology. Our experiments show that the number of iterations required for convergence of CG-QSC (Conjugate Gradient applied to Quadratic Spline Collocation equations) grows linearly with the square root of the number of equations. We implemented the CG and PCG methods for the solution of the Quadratic Spline Collocation (QSC) equations on the iPSC/2 hypercube and present performance evaluation results for up to 32 processors configurations. Our experiments show efficiencies of the order of 90%, for both the fixed and scaled speedups.

Title: A Domain Decomposition Spline Collocation Method for Elliptic Partial Differential Equations
Author(s):C. C. Christara and E. N. Houstis
Date: March 1989
Pages: 7
Note: Proceedings of the fourth Conference on Hypercubes, Concurrent Computers and Applications (HCCA4), March 6-8, 1989, Monterey, CA, U.S.A., pp. 1267-1273.
Abstract:
We consider the parallel solution of a discretization system obtained by approximating the solution of the second order elliptic boundary value problem
Lu == a uxx + b uxy + c uyy + d ux + e uy + fu = g in OMEGA == [ax,bx] X [ay,by] (1.1)
Bu == alpha u + beta un = g_0 on partial OMEGA == boundary of OMEGA (1.2)
using the quadratic spline collocation method. The input functions in the partial differential equation (PDE) problem (1.1), (1.2) are assumed to be functions of $x,y$ in $C^1 [OMEGA]$. In [Hous88], we have studied the computational performance of various direct and iterative methods for solving such discrete systems on sequential machines. An overview of spline collocation methods is presented in Section 2. The objective of this paper is to develop and analyze parallel solvers on shared and non-shared memory MIMD architectures for solving elliptic quadratic spline collocation equations. The class of domain decomposition or substructuring methods has recently been applied by many researchers to solve efficiently PDE discretization systems. Their build in natural parallelism can be easily mapped to various MIMD architectures. Our focus in this study is limited to a parallel formulation of the capacitance method for spline collocation equations. According to the basic idea of this method, the computational effort of solving the original system is reduced to the solution of the so called capacitance matrix problem. For positive definite systems, the capacitance system is solved with preconditioning conjugate gradient (PCG). In the case of quadratic spline collocation equations, in addition to the PCG method, we have developed a scaled Jacobi scheme for its solution, which is formulated and analyzed in Section 3. The parallelization of the capacitance - Jacobi solver, is presented in Section 4 together with its complexity and performance analysis on a hypercube and bus architectures. Finally, the performance on the MIMD machines (NCUBE/7 and SEQUENT) are reported in Section 5.

Title: A Parallel Spline Collocation - Capacitance Method for Elliptic Partial Differential Equations
Author(s): C. C. Christara, E. N. Houstis and J. R. Rice
Date: July 1988
Pages: 10
Note: Proceedings of the 1988 International Conference on Supercomputing (ICS88), July 1988, Saint-Malo, France, pp. 375-384.
Abstract:
We consider the integration of a domain decomposition technique with a new quadratic spline collocation discretization scheme for solving second order elliptic boundary value problems on rectangles. The domain decomposition method is based on the capacitance matrix technique. Due to the limitations of existing methods for solving the corresponding capacitance problem, we develop and analyze iterative methods for its solution. The optimum partitioning and mapping of the underlying computation is studied on hypercube architectures. A numerical realization of this method is presented on NCUBE/7 (128 processors) and its comparative efficiency is measured. The resulting parallel quadratic spline collocation-capacitance method is seen to be efficient in achieving accurate solutions and in using parallel architectures.

Title: Geometry Decomposition based Methods for solving Elliptic Partial Differential Equations
Author(s): C. C. Christara, A. Hadjidimos, E. N. Houstis, J. R. Rice and E. A. Vavalis
Date: May 1988
Pages: 8
Note: Comp. Methods in Flow Analysis, (H. Niki and M. Kawahara, eds), 1988, Univ. of Okayama, Japan, pp. 175-182.
Abstract:

Title: Quadratic Spline Collocation Methods for two point boundary value problems
Author(s): E. N. Houstis, C. C. Christara and J. R. Rice
Date: April 1988
Pages: 18
Note: International Journal for Numerical Methods in Engineering, 26:4, April 1988, pp. 935-952.
Abstract:
A new collocation method based on quadratic splines is presented for second order two point boundary value problems. First, $O(h^4)$ approximations to the first and second derivative of a function are derived using a quadratic spline interpolant of $u$. Then these approximations are used to define an $O(h^4)$ perturbation of the given boundary value problem. Second, the perturbed problem is used to define a collocation approximation at interval midpoints for which an optimal $O(h^{3-j})$ global estimate for the $j$th derivative of the error is derived. Further, $O(h^{4-j})$ error bounds for the $j$th derivative are obtained for certain superconvergence points. It should be observed that standard collocation at midpoints gives $O(h^{2-j})$ bounds. Results from numerical experiments are reported that verify the theoretical behavior of the method.

Title: Performance of Scientific Software
Author(s): E. N. Houstis, J. R. Rice, C. C. Christara and E. A. Vavalis
Date: 1988
Pages: 23
Note: Mathematical Aspects of Scientific Software, The IMA Volumes in Mathematics and its applications, Vol 14, Springer Verlag, 1988, pp. 123-155.
Abstract:
We review performance methodologies used for the evaluation of scientific software in von Neumann architectures. A prototype evaluation facility for second order elliptic partial differential equation (PDE) solvers is described which realizes the main objectives of these methodologies. Finally, the results of an evaluation study for a new class of spline collocation solvers for elliptic PDEs are presented.

Title: Scientific Computing by Numerical Methods
Author(s): C. C. Christara and K. R. Jackson
Date: 1996
Pages: 121
File: survey.pdf
Note: A different version of this article appeared in Encyclopaedia of Applied Physics, Vol 17, 1996, pp. 1-79, and in Mathematical Tools for Physicists.
Abstract:
Numerical methods are an indispensable tool in solving many problems that arise in applied physics. In this article, we briefly survey a few of the most common mathematical problems and review some numerical methods to solve them.
As can be seen from the table of contents, the topics covered in this survey are those that appear in most introductory numerical methods books. However, our discussion of each topic is more brief than is normally the case in such texts and we frequently provide references to more advanced topics that are not usually considered in introductory books.
We originally intended to include a section on the computation of eigenvalues and eigenvectors of matrices, but, because of time and space constraints, we omitted it. The interested reader should consult an introductory numerical methods book, such as those listed in the last section, or an advanced text, such as \cite{Golub-Van-Loan, Parlett, Saad, Stewart, Wilkinson}.

Supervised students' theses

Title: Bilateral XVA pricing under stochastic default intensity: PDE modelling and computation,
Author: Yuwei Chen
Date: February 2023
Pages: 123
Note:Ph.D. thesis, Computer Science Department, University of Toronto
File: ywchen-23-phd.pdf

Abstract:

Adjusting derivative prices to take into account default risk has attracted the attention of several researchers and practitioners, especially after the 2007-2008 financial crisis. The thesis is a study, via numerical Partial Differential Equation (PDE) approaches, of the modeling and computation of valuation adjustment problems in financial derivative pricing if we consider the bilateral counterparty default intensities.
Under some assumptions, the valuation of financial derivatives, including a value adjustment to account for bilateral default risk (XVA), assuming constant default intensities gives rise to a nonlinear PDE \cite{burgard2011partial}. We formulate a penalty iteration method to handle the nonlinearity in this PDE, study its convergence, and extend it to American derivatives, resulting in a double-penalty iteration. We also propose boundary conditions and their discretization for the XVA PDE problem in the cases of call, put options, and forward contracts. Numerical results demonstrate the effectiveness of our methods.
Then, we derive a novel PDE for derivative pricing including XVA, assuming that the default risk of one of the counterparties follows a Cox-Ingersoll-Ross (CIR) process, while the other party has constant default risk. The time-dependent PDE derived is of Black-Scholes type and involves two ``space'' variables, namely the asset price and the buyer default intensity, as well as a nonlinear source term. We formulate boundary conditions appropriate for the default intensity variable. The numerical solution of the PDE uses a penalty-like iteration for handling the nonlinearity. We also develop a novel asymptotic approximation formula for the adjusted price of derivatives, resulting in a very efficient, accurate, and convenient for practitioners formula. Numerical results indicate stable second order convergence for the 2D PDE solution in terms of the discretization size and convergence of order at least 1.5 for the asymptotic approximation in terms of inverse of the mean-reversion rate.
We extend the PDE model we developed to price European derivatives including XVA, considering stochastic counterparty default intensity, to American derivatives. We also extend the asymptotic approximation to the American put option problem. A key step is to derive the asymptotic approximation to the free boundary. We present numerical experiments in order to study the accuracy and effectiveness of the 2D PDE and asymptotic approximations.

Title: Penalty Methods for Nonlinear Problems in Financial Option Pricing
Authors: Ruining (Ray) Wu
Date: February 2021
Pages: 121
Note: MSc Research Paper, Computer Science Department, University of Toronto
File: rwu-21-msc.pdf

Abstract: We study the pricing of nonlinear problems in computational finance by numerical Partial Differential Equation (PDE) methods. We examine the numerical solution obtained via the use of Policy Iteration methods from Hamilton-Jacobi-Bellman equations and introduce new penalty-like iterative methods to solve the nonlinear equations. We consider the following problems: Unequal Borrowing and Lending rates, Stock Borrowing Fees, Stock Borrowing Fees with American exercise rights, Uncertain Volatility, and Transaction cost models. We demonstrate second-order convergence of the solution and of the Greeks, and for the nonlinear problems we study the number of nonlinear iterations taken as a proxy to the total computational cost. Furthermore, the effects of various details such as different boundary conditions and nonuniform discretization grids are studied. Finally, where applicable, we prove monotonicity of the new penalty-like methods that we have introduced and demonstrate that our numerical implementations uphold this property.

Title: Numerical Methods for Pricing Multi-Asset Options
Author: Yuwei Chen
Date: January 2018
Pages: 75
Note:M.Sc. thesis, Computer Science Department, University of Toronto
File: ywchen-18-msc.pdf

Abstract:

We consider the pricing of two-asset European and American options by numerical Partial Differential Equation (PDE) methods, and compare the results with certain analytical formulae. Two cases of options are tested: exchange option and spread option. For exchange options, the analytical formula considered is the (exact) Margrabe formula. For spread options, we consider Kirk's formula and the formula by Li, Deng and Zhou. In pricing European two-asset options, the basic numerical PDE model is the two-dimensional Black-Scholes PDE. Different boundary conditions are considered, and the effect of them on the solution at various points of the grid is studied. Furthermore, various types of non-uniform grids are considered, aiming at reducing the error at certain areas of the grid. Moreover, the effect of the truncated domain on the PDE approximation is studied. We also discuss the effect of certain problem parameters, such as the length of maturity time, and the values of volatility and correlation, to the accuracy and convergence of the PDE approximations. The experiments indicate that the numerical PDE computed price and Greeks are second-order, for appropriately chosen grid discretizations. In the American pricing problem, the discrete penalty method is applied to the linear complementarity problem involving the two-dimensional Black-Scholes PDE and additional constraints. The convergence of the American Spread put option approximation computed with penalty iteration remains second-order, with the number of penalty iterations per timestep remaining small (2-3). We also consider an iterative method with preconditioning techniques for solving the arising large sparse linear system at each timestep, and show that this solution technique is asymptotically optimal.

Title: PDE option pricing: analysis and application to stochastic correlation,
Author: Nat Chun-Ho Leung
Date: June 2017
Pages: 119
Note:Ph.D. Thesis, Depart. of Computer Science, University of Toronto, 2017.
File: nat-17-phd.pdf
Keywords: quantization error, Fourier analysis, digital options, call and put options, stochastic correlation, spread options, quanto options, max options, partial differential equation, localization, boundary conditions, convergence, stability.

Abstract: This thesis is a study of numerical Partial Differential Equation (PDE) methods in financial derivatives pricing. The first part of the thesis is concerned with the behaviour of a numerical PDE solution when the initial condition is not smooth. The second part of the thesis develops computational PDE methods for option pricing problems with stochastic correlation.
In the first part of this thesis, we provide an analysis of the error arising from a non-smooth initial condition when solving a pricing problem with a finite difference method. We build our framework on the sharp error estimate in Giles and Carter (2006), and study three types of non-smoothness that are of financial interest. We show that the error of the numerical solution under Crank-Nicolson-Rannacher timestepping with central spatial differences can be decomposed into two components, respectively a second order error resulting from the approximation to the heat kernel by a discrete operator, and a quantization error that depends on the positioning of non-smoothness on the grid. We discuss how this positioning affects the quality of the numerical solution, and the possibility of an optimal placement of the non-smoothness in the mesh. We also study explicitly the effect of smoothing on the approximation error.
The second part of the thesis focuses on the pricing of European options using a stochastic correlation model. We derive a time-dependent PDE for the pricing problem under stochastic correlation, and develop computational approaches for its solution. The first approach is a finite difference scheme. We study such issues as localization of domain, boundary conditions and stability of the numerical scheme. The second approach is an asymptotic solution of the PDE, appropriate for cases when the correlation process exhibits fast mean reversion and when a numerical PDE solution is considered costly.
Numerical experiments demonstrate the effectiveness of our methods, and the agreement among the two solutions and Monte Carlo simulations. We also experimentally demonstrate the effect of smoothing on the numerical solution, and study the effect of certain problem parameters on the approximate solution.

Title: Efficient and Accurate Numerical PDE Methods for Pricing Financial Derivatives
Author: Mufan (Bill) Li
Date: April 2015
Pages: 35
Note:B.A.Sc. Thesis, Depart. of Electrical and Computer Eng., Univ. of Toronto, 2015.
File: mufanLi-15-basc.pdf
Keywords: American options, (implicit) penalty iteration, operator splitting (explicity penalty) method, convergence, accuracy, complexity, efficiency, Greek calculation.

Abstract: The main difficulty in pricing American options comes from the early exercise right, creating a non-linear constraint on the Black-Scholes PDE. Under a finite difference discretization of the PDE, the price of an American can be approximated, with several techniques to properly handle the American Constraint. While both an iterative penalty method and a direct operator splitting method are convergent, the efficiency and quality require a comprehensive study. Using Crank-Nicolson time stepping and non-uniform grids, the methods are compared in numerical experiments. The criteria include order of convergence, efficiency, and complexity.

Title: Modeling Multi-Factor Financial Derivatives by a Partial Differential Equation Approach with Efficient Implementation on Graphics Processing Units
Author: Duy Minh Dang
Date: August 2011
Pages: 229
Note:Ph.D. Thesis, Computer Science Dept., Univ. of Toronto, 2011.
File: Duy_Minh_Dang_PhD_Thesis.pdf
Keywords: cross-currency interest rate derivatives, Power Reverse Dual Currency (PRDC) swaps, Bermudan cancelability, knockout, and FX Target Redemption, Alternating Direction Implicit (ADI) timestepping methods, multi-GPU platforms, multi-asset American options.

Abstract: This thesis develops efficient modeling frameworks via a Partial Differential Equation (PDE) approach for multi-factor financial derivatives, with emphasis on three-factor models, and studies highly efficient implementations of the numerical methods on novel high-performance computer architectures, with particular focus on Graphics Processing Units (GPUs) and multi-GPU platforms/clusters of GPUs. Two important classes of multi-factor financial instruments are considered: cross-currency/foreign exchange (FX) interest rate derivatives and multi-asset options.

For cross-currency interest rate derivatives, the focus of the thesis is on Power Reverse Dual Currency (PRDC) swaps with three of the most popular exotic features, namely Bermudan cancelability, knockout, and FX Target Redemption. The modeling of PRDC swaps using one-factor Gaussian models for the domestic and foreign interest short rates, and a one-factor skew model for the spot FX rate results in a time-dependent parabolic PDE in three space dimensions. Our proposed PDE pricing framework is based on partitioning the pricing problem into several independent pricing subproblems over each time period of the swap's tenor structure, with possible communication at the end of the time period. Each of these subproblems requires a solution of the model PDE. We then develop a highly efficient GPU-based parallelization of the Alternating Direction Implicit (ADI) timestepping methods for solving the model PDE. To further handle the substantially increased computational requirements due to the exotic features, we extend the pricing procedures to multi-GPU platforms/clusters of GPUs to solve each of these independent subproblems on a separate GPU. Numerical results indicate that the proposed GPU-based parallel numerical methods are highly efficient and provide significant increase in performance over CPU-based methods when pricing PRDC swaps. An analysis of the impact of the FX volatility skew on the price of PRDC swaps is provided.

In the second part of the thesis, we develop efficient pricing algorithms for multi-asset options under the Black-Scholes-Merton framework, with strong emphasis on multi-asset American options. Our proposed pricing approach is built upon a combination of (i) a discrete penalty approach for the linear complementarity problem arising due to the free boundary and (ii) a GPU-based parallel ADI Approximate Factorization technique for the solution of the linear algebraic system arising from each penalty iteration. A timestep size selector implemented efficiently on GPUs is used to further increase the efficiency of the methods. We demonstrate the efficiency and accuracy of the proposed GPU-based parallel numerical methods by pricing American options written on three assets.

Title: A study of a discrete element method based granular dynamics solver
Author: Tina Siu Yee
Date: October 2009
Pages: 43
Note:M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2009.
File: ccc/tyee-09-msc.pdf
Keywords: molecular dynamics, sand, powder, particle simulation, ODE solver, Runge-Kutta, friction, dampening, piling, stiction

Abstract: Granular dynamics is the dynamics of a large set of small particles (grains). Convincing simulation of natural granular phenomena (e.g. dust, sand or powders) is a challenging mathematical and computational problem. Our observation is that the more realistically the collection of grains approaches its static state, the more natural the simulation appears. This study focuses on the simulation of sets of grains as the set approaches its static state. The method begins with a discrete element (also referred to as molecular dynamics) model of the inter-particle contacts within the granular collection. Inertia terms (friction/dampening) are added to the basic contact model to facilitate static piling. An examination of the different contact models on the formation of the final static state and a discussion of the numerical consequences of each model is presented. The discrete element approach demonstrates to be a straightforward and natural way to model many granular behaviors. Its versatility makes it possible to use it to build a general-purpose granular solver.

Title: Multigrid and Spline Collocation Methods on Non-uniform Grids
Author: Guoyu Wang
Date: September 2008
Pages: 98
Note:M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2008.
File: Please contact ccc SYMBOL_AT cs DOT toronto DOT edu
Keywords: quadratic splines, cubic splines, fourth order of convergence, collocation, multigrid, extension operator, restriction operator, preconditioning

Abstract: This thesis describes numerical methods for the solution of one-dimensional second-order differential equations on non-uniform grids. The methods combine the multigrid and spline collocation methodologies. Multigrid methods are presented. Appropriate restriction and extension operators are developed for quadratic and cubic splines on non-uniform grids. Multigrid methods are applied to quadratic and cubic spline collocation equations arising from the discretization of one-dimensional second-order differential equations. The convergence analysis of a two-grid method integrated with a damped Richardson relaxation scheme as smoother shows that the rate of convergence is faster than $\frac{1}{2}$ for a model second-order equation with homogeneous Dirichlet boundary conditions. An adaptive technique for estimating the minimal acceptable coarse grid size is proposed. Numerical experiments using full multigrid methods are performed on various one-dimensional second-order differential equations discretized on non-uniform grids, and the conditions under which they are convergent and efficient are studied. The numerical results confirm our analysis.

Title: Bi-quartic Spline Collocation Methods for Fourth-order Boundary Value Problems with an Application to the Biharmonic Dirichlet Problem
Author: Jingrui Zhang
Date: January 2008
Pages: 160
Note:Ph.D. Thesis, Computer Science Dept., Univ. of Toronto, 2008.
File: ccc/jingrui-08-phd.ps.gz
Keywords: sixth order of convergence, quartic splines, collocation, fast Fourier Transform (FFT), GMRES, preconditioning

Abstract: We present optimal bi-quartic spline collocation methods for the numerical solution of fourth-order boundary value problems on rectangular domains, and apply a particular instance of these methods to the biharmonic Dirichlet problem.

The methods are based on quartic splines and the collocation discretization methodology with midpoints of a uniform partition as the collocation points. While the standard quartic spline method usually provides second-order approximations, the two bi-quartic spline collocation methods developed in this thesis are observed to produce approximations which are of sixth order at grid points and midpoints, and of fifth order at other points. Both are based on high order perturbations of the differential and boundary conditions operators. The one-step (extrapolated) method forces the approximations to satisfy a perturbed problem, while the three-step (deferred-correction) method proceeds in three steps and perturbs the right sides of the linear equations appropriately.

The three-step bi-quartic spline collocation method is then applied to the biharmonic Dirichlet problem and a fast solver for the resulting linear systems is developed. The fast solver consists of solving an auxiliary biharmonic problem with Dirichlet and second derivative boundary conditions along the two opposite boundaries, and a one-dimensional problem related to the two opposite boundaries. The linear system arising from the bi-quartic spline collocation discretization of the auxiliary problem is solvable by fast Fourier transforms. The problem related to the two opposite boundaries is solved by preconditioned GMRES (PGMRES). We present a detailed analysis of the performance of PGMRES by studying bounds for the eigenvalues of the preconditioned matrix. We show that the number of PGMRES iterations is independent of the grid size $n$. As a result, the complexity of the fast solver is $O(n^2\log(n))$. Numerical experiments from a variety of problems, including practical applications and problems more general than the analysis assumes, verify the accuracy of the discretization scheme and the effectiveness of the fast solver.

Title: Adaptive Finite Difference Methods for Valuing American Options
Author: Duy Minh Dang
Date: September 2007
Pages: 128
Note:M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2007.
File: dmdang-07-msc.pdf
Keywords: adaptive mesh selection, monitor function, finite differences, Crank-Nicolson, European option, American option, penalty method, projected SOR

Abstract: We develop space-time adaptive methods for valuing American options with strong emphasis on American put options. We examine the application of adaptive techniques to the Black-Scholes partial differential equation problem associated with an American put option in the context of non-uniform second-order finite differences. At certain timesteps, we obtain a redistribution of the spatial points based on a monitor function that attempts to equidistribute the error. The proposed finite difference discretization on non-uniform grids and redistribution of the spatial points lead to linear complementarity problems with $\Mcal$-matrices. The Projected Successive Over-relaxation and a penalty method are considered to handle the free boundaries. We study and compare the accuracy and efficiency of the considered methods. A complete proof of convergence and uniqueness of the projected SOR method under certain conditions is also presented.

Title: Quartic Spline Collocation Methods For Second-Order Two-Point Boundary Value ODE Problems
Author: Guohong Liu
Date: August 2007
Pages: 130
Note:M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2007.
File: ghliu-07-msc.ps.gz
Keywords: sixth order of convergence, quartic splines, collocation, boundary conditions

Abstract: Collocation methods based on quartic splines are presented for second-order two-point boundary value problems. In addition to the boundary conditions specified by the problem, extra boundary conditions are introduced in order to uniquely determine the quartic spline approximation. The standard quartic spline collocation method gives fourth order bounds. Two optimal methods, namely the extrapolated (one-step) and the deferred-correction (two-step) methods, are formulated based on appropriate extra boundary conditions and an appropriate perturbation of the operators of the differential equation, boundary conditions and extra boundary conditions. The convergence analysis and error bounds are developed using a Green's functions approach. The analysis shows that the maximum discrete error is of sixth order, and the maximum global error is of fifth order for the optimal methods.
The properties of the matrices arising from the optimal methods for a certain class of problems are studied. Non-optimal collocation methods based on different extra boundary conditions are also investigated. Problems with layers are also considered, and a grid refinement technique is presented. The theoretical behavior of the collocation methods is verified by numerical results.

Title: Pricing Convertible Bonds with Dividend Protection subject to Credit Risk Using a Numerical PDE Approach
Author: Qingkai Mo
Date: April 2006
Pages: 92
Note:M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2006.
File:ccc/moqk-06-msc.pdf
Keywords: convertible bond, dividend protection, credit risk, Conversion Ratio Adjustment, Dividend Pass-Thru, finite differences, Crank-Nicolson, Projected Successive Over-Relaxation (PSOR), penalty method.

Abstract: We develop a pricing model for convertible bonds with dividend protection subject to credit risk by extending the models developed by Tsiveriotis and Fernandes (TF), and by Ayache, Forsyth and Vetzal (AFV). We consider two techniques to incorporate the dividend protection feature: Conversion Ratio Adjustment and Dividend Pass-Thru. We apply finite difference methods to discretize the PDEs associated with our models, and study the Projected Successive Over-Relaxation (PSOR) and penalty methods to handle the free boundaries. We compare these two methods in terms of convergence rate, number of iterations per time step and computation time for pricing convertible bonds without dividends. Finally, we apply the penalty method, the better of the two methods, to solve the systems arising from our models for convertible bonds with dividend protection. We examine the convergence rates and discuss the difference between the results from the extended TF and AFV models, with both dividend protection techniques.

Title: An Efficient Algorithm Based on Quadratic Spline Collocation and Finite Difference Methods for Parabolic Partial Differential Equations
Author: Tong Chen
Date: September 2005
Pages: 92
Note:M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2005.
File:ccc/tchen-05-msc.pdf
Keywords: semi-explicit semi-implicit method, optimal order of convergence, Crank-Nicolson, relaxed conditions for stability

Abstract: An efficient algorithm which combines quadratic spline collocation methods (QSC) for the space discretization and classical finite difference methods (FDMs), such as Crank-Nicolson, for the time discretization to solve general linear parabolic partial differential equations has been studied. By combining QSC and finite differences, a form of the approximate solution of the problem at each time step can be obtained; thus the value of the approximate solution and its derivatives can be easily evaluated at any point of the space domain for each time step.
There are two typical ways for solving this problem: (a) using QSC in its standard formulation, which has low accuracy $\mathcal{O}(h^2)$ and low computational work. More precisely, it requires the solution of a tridiagonal linear system at each time step; (b) using optimal QSC, which has high accuracy $\mathcal{O}(h^4)$ and requires the solution of either two tridiagonal linear systems or an almost pentadiagonal linear system at each time step. A new technique is introduced here which has the advantages of the above two techniques; more precisely, it has high accuracy $\mathcal{O}(h^4)$ and almost the same low computational work as the standard QSC.

Title: Pricing Convertible Bonds using Partial Differential Equations
Author: Lucy Xingwen Li
Date: September 2005
Pages: 81
Note:M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2005.
File:lucy-05-msc.pdf
Keywords: Coupon payments, call, put, convert, discontinuities, implicit methods, successive over-relaxation, penalty method

Abstract:

A Convertible Bond (CB) is a corporate debt security that gives the holder the right to exchange future coupon payments and principal repayment for a prescribed number of shares of equity. Thus, it has both an equity part and a fixed-income part, and may contain some additional features, such as callability and puttability.
In this paper, we study the model for valuing Convertible Bonds with credit risk originally developed by Kostas Tsiveriotis and Chris Fernandes (TF). The Convertible Bond is a derivative of the stock price, and the pricing model developed by TF is based on a free boundary value problem associated with a pair of parabolic Partial Differential Equations (PDEs) with discontinuities at the time points when there is a coupon payment, or when the bond is converted, or when it is called back (purchased) by the issuer, or when it is put (sold) to the issuer. We explore the possible derivation of the TF model and study the convergence of several numerical methods for solving the free boundary value problem associated with the TF model. In particular, we consider the Successive Over-Relaxation (SOR) iteration and a penalty method for solving the linear complementarity problem used to handle the free boundary. Special emphasis is given to the effectiveness of the numerical scheme.

Title: Spline Collocation on Adaptive Grids and Non-Rectangular Domains
Author:Kit Sun Ng
Date: June 2005
Pages: 187
Note:Ph.D. Thesis, Computer Science Dept., Univ. of Toronto, 2005.
File:ngkit-05-phd.ps.gz
Keywords: quadratic splines, cubic splines, optimal convergence, non-uniform grid, adaptive methods, gridsize estimator, error estimator, L-shaped domains, T-shaped domains, skipped grid, moving mesh PDE.

Abstract:

Optimal spline collocation is a computationally cost effective finite element method that has been relatively recently developed. Although there are many important results for optimal spline collocation, we believe the method can further be extended to solve a series of more relevant and difficult problems.

We first consider using optimal spline collocation on non-uniform partitions to solve one-dimensional linear second-order Boundary Value Problems (BVPs) that have solutions with regions of high variations, such as boundary layers. Optimal quadratic spline collocation (QSC) is extended to non-uniform partitions by using a mapping function from uniform to non-uniform partitions. Optimal cubic spline collocation (CSC) on non-uniform partitions is also examined. The formulation and analysis of the CSC methods are based, as in the QSC case, on a mapping function from uniform to non-uniform partitions. However, this method is converted into a mapping-free method and is integrated with adaptive techniques that include error and grid size estimators.

Spline collocation is then considered for solving two-dimensional linear second-order BVPs. We extend optimal spline collocation to solving BVPs on L- and T-shaped domains, which occur often in practice. The CSC method can handle general boundary conditions, which had previously never been done, but our approach requires a third step to achieve optimal convergence. We also consider two adaptive approaches for solving problems on rectangular domains whose solutions have rapid variations. As in the one-dimensional case, such problems present a significant challenge to obtain accurate approximations at minimal computational cost. We first implement a static method, which involves adding collocation points to regions where the solution has sharp variations, with CSC that uses skipped grids, which are grids that have several levels of refinement, in certain regions of the domain. An optimal skipped grid method is developed by converting the second-order BVP into a system of three first-order BVPs. We then implement a dynamic method, which involves moving the collocation points to regions where the solution has sharp variations, with CSC that incorporates the {\it Moving Mesh Partial Differential Equation} (MMPDE). The MMPDE method is further integrated into an adaptive method that can solve BVPs within a specific error tolerance.

Title: Multigrid and Cubic Spline Collocation Methods for Advection Equations
Author: Zheng Zeng
Date: April 2005
Pages: 111
Note:M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2005.
File: ccc/zengz-05-msc.ps.gz
Keywords: Shallow water equations, semi-Lagrangian methods, algebraic elimination, analytic elimination, systems of PDEs

Abstract:
This thesis describes numerical methods for the solution of advection equations in one-dimensional space. The methods are based on combining the multigrid and cubic spline collocation (CSC) methodologies.
Multigrid methods for cubic splines are presented. Appropriate restriction and extension operators are developed for cubic splines that satisfy various boundary conditions (BCs), i.e., Dirichlet, Neumann, periodic and general BCs. The multigrid methods are applied to CSC equations arising from the discretization of one-dimensional second-order differential equations. The convergence analysis of a two-grid method integrated with a damped Richardson relaxation scheme as smoother shows the rate of convergence is equal to or faster than $1/2$ for a model second order equation with periodic BCs.
Multigrid methods for cubic splines are then extended to the solution of one-dimensional shallow water equations (SWEs). The SWEs are discretized in time with a semi-Lagrangian semi-implicit scheme and in space with CSC methods. At each time step, there are three different space discretization approaches: Analytic Elimination - Discretization (AED), Discretization -Algebraic Elimination (DAE), and Discretization - No Elimination (DNE). We discuss advantages and disadvantages of each, and develop new numerical methods based on multigrid and CSC for solving the SWEs discretized with either the AED or the DNE approaches. We compare our methods with Jacobi's iterative method applied to the same problems. Through comparison and convergence discussion, we show that our multigrid methods for CSC are convergent and efficient. The numerical results confirm our analysis.

Title: Towards a high-performance method for the biharmonic Dirichlet problem
Authors: Jingrui Zhang
Date: April 2004
Pages: 77
Note:Research paper (equivalent to M.Sc. Thesis), Computer Science Dept., Univ. of Toronto, 2004.
File:ccc/jingrui-04-msc.ps.gz
Keywords: spline collocation; fourth-order two-point boundary value problem; optimal order of convergence; FFT; GMRES; preconditioning; eigenvalues; eigenvectors

Abstract:
The biharmonic Dirichlet problem appears in many applications. Therefore, there is a lot of interest in developing high-performance methods for its solution. Two important components of performance are accuracy and efficiency. In this research, we combine a sixth order discretization method and a Fast Fourier Transform (FFT)-based solver for the solution of the biharmonic Dirichlet problem.

The discretization of the problem uses the collocation methodology and quartic splines. Quartic Spline Collocation (Quartic SC) methods that produce sixth order approximations have been recently developed for the solution of fourth order Boundary Value Problems (BVPs). In this paper, we apply an adjustment to the quartic spline basis functions so that the Quartic SC matrix has more favorable properties. Particular attention is given to the one-dimensional harmonic problem. We study the properties of two Quartic SC matrices, more specifically, one that corresponds to Dirichlet and Neumann conditions (the standard harmonic problem) and another that corresponds to Dirichlet and second derivative conditions. The latter is solvable by FFTs and used as preconditioner for the first. We develop formulae for the eigenvalues and eigenvectors of the preconditioned matrix. We show that three iterations suffice to obtain convergence for the GMRES preconditioned method. Numerical experiments verify the effectiveness of the solver and preconditioner. We also discuss the application of the preconditioned GMRES method to the two-dimensional biharmonic Dirichlet problem.

Title: Quartic-Spline Collocation Methods for Fourth-Order Two-Point Boundary Value Problems
Author:Ying Zhu
Date: April 2001
Pages: 73
Note: M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2001.
File: ccc/yz-01-msc.ps.gz

Abstract:
This thesis presents numerical methods for the solution of general linear fourth-order boundary value problems in one dimension. The methods are based on quartic splines, that is, piecewise quartic polynomials with C^3 continuity, and the collocation discretization methodology with the midpoints of a uniform partition being the collocation points.
The standard quartic-spline collocation method is shown to be second order. Two sixth-order quartic-spline collocation methods are developed and analyzed. They are both based on a high order perturbation of the differential equation and boundary conditions operators. The one-step method forces the approximation to satisfy a perturbed problem, while the three-step method proceeds in three steps and perturbs the right sides of the equations appropriately.
The error analysis follows the Green's function approach and shows that both methods exhibit optimal order of convergence, that is, they are locally sixth order on the gridpoints and midpoints of the uniform partition, and fifth order globally. The properties of the matrices arising from a restricted class of problems are studied. Analytic formulae for the eigenvalues and eigenvectors are developed, and related to those arising from quadratic-spline collocation matrices. Numerical results verify the orders of convergence predicted by the analysis.

Title: High-Order Spatial Discretization Methods for the Shallow Water Equations
Author:Anita W. Tam
Date: February 2001
Pages: 160
Note: Ph.D. Thesis, Computer Science Dept., Univ. of Toronto, 2001.
File: tam-01-phd.ps.gz

Abstract:
We present new numerical methods for the shallow water equations on a sphere in spherical coordinates. In our implementation, the equations are discretized in time with the two-level semi-Lagrangian semi-implicit (SLSI) method, and in space on a staggered grid with the quadratic spline Galerkin (QSG) and the optimal quadratic spline collocation (OQSC) methods. When discretized on a uniform spatial grid, the solutions are shown through numerical experiments to be fourth-order in space locally at the nodes and midpoints of the spatial grids, and third-order globally.

We also show that, when applied to a simplified version of the shallow water equations, each of our algorithms yields a neutrally stable solution for the meteorologically significant Rossby waves. Moreover, we demonstrate that the Helmholtz equation associated with the shallow water equations should be derived algebraically rather than analytically in order for the algorithms to be stable with respect to the Rossby waves. These results are verified numerically using Boyd's equatorial wave equations with initial conditions to generate a soliton.

We then analyze the performance of our methods on various staggered grids --- the A-, B-, and C-grids. From an eigenvalue analysis of our simplified version of the shallow water equations, we conclude that, when discretized on the C-grid, our algorithms faithfully capture the group velocity of inertia-gravity waves. Our analysis suggests that neither the A- nor B-grids will produce such good results. Our theoretical results are supported by numerical experiments, in which we discretize Boyd's equatorial wave equations using different staggered grids and set the initial conditions for the problem to generate gravitation modes instead of a soliton. With both the A- and B-grids, some waves are observed to travel in the wrong direction, whereas, with the C-grid, gravity waves of all wavelengths propagate in the correct direction.

Title: Quadratic Spline Collocation Methods for Systems of Elliptic PDEs
Authors:Kit Sun Ng
Date: April 2000
Pages: 67
Note: M.Sc. Thesis, Computer Science Departmnent, University of Toronto, 2000.
File: ccc/ngkit-00-msc.ps.gz
Abstract:
We consider Quadratic Spline Collocation (QSC) methods for solving systems of two linear second-order PDEs in two dimensions. Optimal order approximation to the solution is obtained, in the sense that the convergence order of the QSC approximation is the same as the order of the quadratic spline interpolant.

We study the matrix properties of the linear system arising from the discretization of systems of two PDEs by QSC. We give sufficient conditions under which the QSC linear system is uniquely solvable and the optimal order of convergence for the QSC approximation is guaranteed. We develop fast direct solvers based on Fast Fourier Transforms (FFTs) and iterative methods using multigrid or FFT preconditioners for solving the above linear system.

Numerical results demonstrate that the QSC methods are fourth order locally on certain points and third order globally, and that the computational complexity of the linear solvers developed is almost asymptotically optimal. The QSC methods are compared to conventional second order discretization methods and are shown to produce smaller approximation errors in the same computation time, while they achieve the same accuracy in less time.

Title: Numerical Solution of the Shallow-Water Equations on Distributed Memory Systems
Author:Xiaoliang (Lloyd) Ding
Date: October 1998
Pages: 73
Note: M.Sc. Thesis, Computer Science Departmnent, University of Toronto, 1998.
File: ding.98.msc.ps.gz
Abstract:
The shallow-water equations are often used as a mathematical model when numerical methods for solving weather or climate prediction problems are tested. This thesis studies the performance and scalability of numerical methods for the shallow-water equations on distributed memory systems. Time integration of the numerical methods is based on a two-time-level semi-implicit semi-Lagrangian scheme. To solve the Helmholtz problem that arises at each time-step, a fast direct solver based on FFTs is used. This technique requires a data transposition, which is a time consuming operation on a distributed memory system.

The data transposition requires an all-to-all type communication. The Direct Exchange algorithm, commonly used for this task, is inefficient if the number of processors is large. We investigate a series of more sophisticated techniques: the Ring Exchange, Mesh Exchange and Cube Exchange algorithms. We compare the performance of these techniques with that of the Message Passing Interface collective communication routine. We perform an isoefficiency analysis to measure the scalability of the parallel system. The numerical experiments are carried out on a Cray T3E with 512 processors and on an ethernet-connected cluster of 36 workstations.

Both the analysis and our numerical results indicate that the Cube Exchange and Mesh Exchange algorithms perform better than the other methods. However, even with the more sophisticated techniques, the parallel algorithm for the shallow-water equations remains essentially unscalable.

Title: Fast Fourier transform solvers for quadratic spline collocation
Author: A. Constas
Date: July 1996
Pages: 64
Note: M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 1996.
File: ccc/constas-96-msc.ps.gz
Abstract:
The solution of linear systems arising from the discretisation of second order elliptic Boundary Value Problems (BVPs) defined in rectangular domains by solvers employing Fast Fourier Transforms (FFTs) is studied. The types of boundary conditions considered are homogeneous Dirichlet, periodic, homogeneous Neumann or combinations of these. Two-dimensional and three-dimensional problems are considered. Quadratic spline collocation methods are used for the discretisation of the BVPs. The FFT solvers are first developed for differential operators with constant coefficients and even order derivatives. Diagonal scalings of these solvers are then used as preconditioners for the solution of linear systems arising from the discretisation of general second order linear elliptic BVPs. Several acceleration methods are tested. Results from numerical experiments that demonstrate the asymptotically optimal convergence and computational performance of the FFT solvers are presented.