Selected publications by C. C. Christara

Supervised students' theses

Title:Adaptive and high-order methods for valuing American options
Authors:Christina 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: Accepted for publication in the Journal of Computational Finance
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 PDE Pricing Framework for Cross-Currency Interest Rate Derivatives
Author:Duy Minh Dang, Christina 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: 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: Please contact ccc SYMBOL_AT cs DOT toronto DOT edu

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:To appear 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: Please contact ccc SYMBOL_AT cs DOT toronto DOT edu

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: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: 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.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: Accepted for publication in the 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: multig.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: 79
Note: Encyclopaedia of Applied Physics, Vol 17, 1996, pp. 1-79.
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}.

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: Please contact ccc SYMBOL_AT cs DOT toronto DOT edu
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: Please contact ccc SYMBOL_AT cs DOT toronto DOT edu
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: Please contact ccc SYMBOL_AT cs DOT toronto DOT edu
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.