Selected publications by Christina Christara

A highorder deferred correction method for the solution of
free boundary problems using penalty iteration, with an
application to American option pricing,
Dawei Wang, Kirill Serkh, and Christina Christara

Bilateral XVA pricing under stochastic default intensity:
PDE modelling and computation,
Yuwei Chen and Christina Christara

Penalty and PenaltyLike Methods for Nonlinear HJB PDEs,
Christina Christara and Ruining Wu

Penalty Methods for Bilateral XVA Pricing,
in European and American Contingent Claims by a PDE Model,
Yuwei Chen and Christina Christara

Analysis of quantization error in financial pricing
via finite difference methods,
Christina C. Christara and Nat ChunHo Leung

Partial Differential Equation pricing of contingent claims
under stochastic correlation,
Nat ChunHo Leung, Christina C. Christara and DuyMinh Dang

Spread option pricing using ADI methods,
Vida HeidarpourDehkordi and Christina C. Christara

PDE option pricing with variable correlations,
Christina C. Christara and Nat ChunHo Leung

Option pricing in jump diffusion models with
quadratic spline collocation,
Christina C. Christara and Nat ChunHo Leung

An efficient numerical PDE approach for pricing foreign exchange
interest rate hybrid derivatives,
Duy Minh Dang, Christina C. Christara, Kenneth R. Jackson and Asif Lakhany

Graphics processing unit pricing
of exotic crosscurrency interest rate derivatives with a
foreign exchange volatility skew model,
Duy Minh Dang, Christina C. Christara and Kenneth R. Jackson

A Highly Efficient Implementation on GPU Clusters of
PDEBased Pricing Methods for
PathDependent Foreign Exchange Interest Rate Derivatives,
Duy Minh Dang, Christina C. Christara and Kenneth R. Jackson

An Efficient GPUBased Parallel Algorithm for
Pricing MultiAsset American Options,
Duy Minh Dang, Christina C. Christara and Kenneth R. Jackson

Adaptive and highorder methods for valuing American options,
Christina C. Christara and Duy Minh Dang

A Parallel Implementation on GPUs of ADI Finite Difference Methods
for Parabolic PDEs with Applications in Finance,
Duy Minh Dang, Christina C. Christara and Kenneth R. Jackson

A PDE Pricing Framework for CrossCurrency Interest Rate Derivatives,
Duy Minh Dang, Christina C. Christara, Kenneth R. Jackson and Asif Lakhany

Pricing MultiAsset American Options on Graphics
Processing Units Using a PDE Approach,
Duy Minh Dang, Christina C. Christara and Kenneth R. Jackson

Quartic spline collocation for secondorder boundary value problems,
Christina C. Christara and Guohong Liu

Quadratic spline collocation for onedimensional
linear parabolic partial differential equations,
Christina C. Christara, Tong Chen and Duy Minh Dang

Quartic spline collocation for fourthorder boundary value problems,
Christina C. Christara, Ying Zhu and Jingrui Zhang

A HighPerformance Method for the Biharmonic Dirichlet Problem on Rectangles,
Christina C. Christara and Jingrui Zhang

Optimal Quadratic and Cubic Spline Collocation on Nonuniform Partitions,
Christina C. Christara, and Kit Sun Ng

Adaptive Techniques for Spline Collocation,
Christina C. Christara, and Kit Sun Ng

GIAinduced secular variations in the Earth's long wavelength
gravity field: Influence of 3D viscosity variations
(short communication),
Konstantin Latychev, Jerry X. Mitrovica, Mark E. Tamisiea, Jeroen Tromp,
Christina C. Christara and Robert Moucha

Glacial isostatic adjustment on 3D Earth models: a finitevolume formulation,
Konstantin Latychev, Jerry X. Mitrovica, Jeroen Tromp, Mark E. Tamisiea,
Dimitri Komatitsch, and Christina C. Christara

Optimal Quadratic Spline Collocation on NonUniform Partitions,
Christina C. Christara, and Kit Sun Ng

Optimal Cubic Spline Collocation on NonUniform Partitions,
Christina C. Christara, and Kit Sun Ng

Quadratic Spline Galerkin Method for the Shallow Water Equations,
Anita T. Layton, Christina C. Christara, and Kenneth R. Jackson

Optimal Quadratic Spline Collocation Methods for the
Shallow Water Equations,
Anita T. Layton, Christina C. Christara, and Kenneth R. Jackson

Quadratic Spline Collocation Revisited:
Extension to Systems of Elliptic PDEs,
Christina C. Christara and Kit Sun Ng

Fast Fourier Transform Solvers and Preconditioners
for Quadratic Spline Collocation,
Christina C. Christara and Kit Sun Ng

An Efficient Transposition Algorithm for Distributed Memory Computers,
C. Christara, X. Ding and K. R. Jackson

Multigrid and Multilevel Methods for Quadratic Spline Collocation,
Christina C. Christara and Barry Smith

Parallel solvers for spline collocation equations, C. C. Christara

Quadratic spline collocation methods for
elliptic partial differential equations, C. C. Christara

High performance computing of elliptic partial differential equations
with spline collocation, C. C. Christara

Parallel computation of partial differential equations
on distributed memory machines, C. C. Christara

Multicolour orderings and iterative methods for elliptic equations,
C. C. Christara

Domain decomposition and incomplete factorisation methods for
partial differential equations, C. C. Christara

Schur complement preconditioned conjugate gradient methods for
spline collocation equations, C. C. Christara

Conjugate gradient methods for spline collocation equations, C. C. Christara

A Domain Decomposition Spline Collocation Method for Elliptic Partial
Differential Equations,
C. C. Christara and E. N. Houstis

A Parallel Spline Collocation  Capacitance Method for Elliptic Partial
Differential Equations,
C. C. Christara, E. N. Houstis and J. R. Rice

Geometry Decomposition based Methods for solving Elliptic Partial
Differential Equations,
C. C. Christara, A. Hadjidimos, E. N. Houstis, J. R. Rice and E. A. Vavalis

Quadratic Spline Collocation Methods for two point boundary value problems,
E. N. Houstis, C. C. Christara and J. R. Rice

Performance of Scientific Software,
E. N. Houstis, J. R. Rice, C. C. Christara and E. A. Vavalis

Scientific Computing by Numerical Methods,
C. C. Christara and K. R. Jackson
Supervised students' theses

Highorder numerical methods for parabolic PDEs
with nonsmooth data: a perspective from option pricing,
Dawei Wang, Ph.D. thesis, University of Toronto,
June 2024, 124 pages.

Numerical PDE Methods for Pricing Bilateral XVA,
Yuwei Chen, Ph.D. thesis, University of Toronto,
February 2023, 123 pages.

Penalty Methods for Nonlinear Problems in Financial Option Pricing,
Ruining (Ray) Wu, MSc Research Paper,
February 2021, 121 pages.

Numerical Methods for Pricing MultiAsset Options,
Yuwei Chen, M.Sc. thesis, University of Toronto,
January 2018, 75 pages.

PDE option pricing: analysis and application
to stochastic correlation,
Nat ChunHo Leung, Ph.D. thesis, University of Toronto,
June 2017, 119 pages.

Efficient and accurate numerical PDE methods
for pricing financial derivatives,
Mufan (Bill) Li, B.A.Sc. thesis, University of Toronto,
April 2015, 35 pages.

Modeling multifactor financial derivatives
by a Partial Differential Equation approach
with efficient implementation on Graphics Processing Units,
Duy Minh Dang, Ph.D. thesis, University of Toronto,
August 2011, 228 pages.

A study of a discrete element method based granular dynamics solver,
Tina Yee, M.Sc. thesis, University of Toronto, 2009, 43 pages.

Multigrid and Spline Collocation Methods on Nonuniform Grids,
Guoyu Wang, M.Sc. thesis, University of Toronto, September 2008, 98 pages.

Biquartic Spline Collocation Methods for Fourthorder
Boundary Value Problems with an Application to
the Biharmonic Dirichlet Problem,
Jingrui Zhang, Ph.D. thesis,
University of Toronto, January 2008, 160 pages.

Adaptive Finite Difference Methods for Valuing American Options,
Duy Minh Dang, M.Sc. Thesis

Quartic Spline Collocation Methods For SecondOrder TwoPoint
Boundary Value ODE Problems,
Guohong Liu, M.Sc. thesis

Pricing Convertible Bonds with Dividend Protection
subject to Credit Risk Using a Numerical PDE Approach,
Qingkai Mo, M.Sc. Thesis

An Efficient Algorithm Based on Quadratic Spline Collocation and
Finite Difference Methods for Parabolic Partial Differential Equations,
Tong Chen, M.Sc. Thesis

Pricing Convertible Bonds using Partial Differential Equations,
Lucy Xingwen Li, M.Sc. Thesis

Spline Collocation on Adaptive Grids and NonRectangular Domains,
Kit Sun Ng, Ph.D. Thesis

Multigrid and Cubic Spline Collocation Methods for Advection Equations,
Zheng Zeng, M.Sc. Thesis

Towards a highperformance method for the biharmonic Dirichlet problem,
Jingrui Zhang, Research Paper (equivalent to M.Sc. Thesis)

QuarticSpline Collocation Methods for FourthOrder TwoPoint
Boundary Value Problems,
Ying Zhu, M.Sc. Thesis

HighOrder Spatial Discretization Methods for the Shallow Water
Equations,
Anita W. Tam, Ph.D. Thesis, cosupervised by Ken Jackson

Quadratic Spline Collocation Methods for Systems of Elliptic PDEs,
Kit Sun Ng, M.Sc. Thesis

Numerical Solution of the ShallowWater Equations on Distributed
Memory Systems,
Xiaoliang (Lloyd) Ding, M.Sc. Thesis, cosupervised by Ken Jackson

Fast Fourier transform solvers for quadratic spline collocation,
Athanasia Constas, M.Sc. Thesis
Title: 
A highorder 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:  A later version published in
Journal of Computational and Applied Mathematics,
Volume 432, November 2023, 115272,
Link to journal page

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 highorder deferred correction algorithm combined with
penalty iteration for solving free and moving boundary problems, using a
fourthorder 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
fourthorder 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 righthand 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
fourthorder 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 highorder 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 236259, March 2023.
link to journal page

Keywords:  Partial Differential Equations, BlackScholes,
CrankNicolson 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 20072008 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
CoxIngersollRoss (CIR) process, while the other party
has constant default risk.
The timedependent PDE derived is of BlackScholes 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 penaltylike 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 PenaltyLike 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 119, 2022.
link to journal page

Keywords:  Partial Differential Equations, BlackScholes,
nonlinear iteration, finite differences, CrankNicolson,
control problem, HamiltonJacobiBellman (HJB) equation,
transaction costs, stock borrowing fees, penalty methods, policy iteration.

File: 
hjbpenalty.pdf

Abstract:
There are numerous financial problems that can be posed as
optimal control problems, leading to HamiltonJacobiBellman
or HamiltonJacobiBellmanIssacs 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 penaltylike method for the use with European exercise rights,
and extend this to American exercise rights resulting
in a doublepenalty method.
We also use our findings to improve the policy iteration algorithms
described in [Forsyth, Labahn 2007].
Numerical results are provided showing clear secondorder 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, BlackScholes,
credit risk, default risk, credit valuation adjustment,
call option, put option, forward contract,
nonlinear iteration, finite differences, CrankNicolson,
control problem, HamiltonJacobiBellman (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 20072008 financial crisis.
Under some assumptions, the valuation of financial derivatives,
including a value adjustment to account for default risk
(the socalled 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 doublepenalty 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 ChunHo Leung 
Date:  July 2017, rev. March 2018

Pages:  27

Keywords: 
nonsmooth initial conditions, option pricing,
numerical solution, partial differential equation,
convectiondiffusion equations, Fourier analysis, finite difference methods

Note:  A version of this paper is published in SIAM J. Numerical Analysis, 56:3, 2018, pp. 17311757

File:  ccc/papererror.pdf

Abstract:
In this paper, we study the error of a second order finite
difference scheme for the onedimensional convectiondiffusion
equation. We consider nonsmooth 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
nonsmoothness on the numerical grid. Based on our analysis,
we study the issue of optimal placement of such
nonsmoothness 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 ChunHo Leung, Christina C. Christara and DuyMinh 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 B1B31 (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 multidimensional timedependent 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 HeidarpourDehkordi and Christina C. Christara 
Date:  July 2017

Pages:  17

Keywords: 
Modified CraigSneyd, Alternating Direction Implicit method,
twodimensional BlackScholes, 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. 353369.

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 twodimensional BlackScholes
Partial Differential Equation (PDE),
use finite difference discretization in space
and consider CrankNicolson (CN)
and Modified CraigSneyd (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 onedimensional 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 higherdimensional problems.
The results from spread option pricing are compared with those
obtained from the closedform approximation formulae
of Kirk (1995), Venkatramanan and Alexander (2011),
Monte Carlo simulations, and
the BrennanSchwartz ADI DouglasRachford method,
as implemented in MATLAB.
In all spread option test cases we considered, including American ones,
our ADIMCS method, implemented on appropriate nonuniform 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 ChunHo 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. 4357.

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 regimeswitching 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 ChunHo Leung 
Date:  May 2014; revised May 2015

Pages:  15

Keywords: 
quadratic spline, collocation, finite difference,
European option, American option, partial integrodifferential equation,
Merton model, Kou model, calculation of Greeks

Note:  Published in
Applied Mathematics and Computation, Vol 279, 10 April 2016, pp 2842.

File:  ccc/paperjump.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: 
PowerReverse DualCurrency (PRDC) swaps, Target Redemption (TARN),
knockout, Partial Differential Equation (PDE), finite differences,
nonuniform grids.

Note:  Published in
Journal of Computational Finance, Vol 18, Issue 4, 2015, pp 3993.

File:  prdcTARN.pdf,
Link to the paper on the SSRN

Abstract:
We discuss efficient pricing methods via a Partial Differential Equation (PDE)
approach for longdated foreign exchange (FX) interest rate hybrids
under a threefactor multicurrency pricing model with FX volatility skew.
The emphasis of the paper is on
PowerReverse DualCurrency (PRDC) swaps with popular exotic features,
namely knockout and FX Target Redemption (FXTARN).
Challenges in pricing these derivatives via a PDE approach
arise from the highdimensionality of the model PDE,
as well as from the complexities in handling the exotic features,
especially in the case of the FXTARN provision,
due to its pathdependency.
Our proposed PDE pricing framework
for FXTARN PRDC swaps 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 pricing subproblems can be viewed as
equivalent to a knockout PRDC swap with a known timedependent barrier,
and requires a solution of the model PDE, which, in our case,
is a timedependent parabolic PDE in three space dimensions.
Finite difference schemes on nonuniform 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 crosscurrency 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.~16091625.

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
crosscurrency interest rate derivatives via a Partial Differential
Equation (PDE) approach. In particular, we focus on the GPUbased
parallel pricing of longdated foreign exchange (FX) interest rate hybrids,
namely Power Reverse Dual Currency (PRDC) swaps
with Bermudan cancelable features.
We consider a threefactor pricing model with
FX volatility skew which results in a
timedependent 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
PDEBased Pricing Methods for
PathDependent 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.~107126,
Proceedings of
The 13th International Conference on Computational Science
and Its Applications (ICCSA 2013),
Ho Chi Minh City, Vietnam, 2427 June 2013.

File:  dangccckrjiccsa13.pdf,
Link to the paper on the SSRN

Abstract:
We present a highly efficient parallelization of the computation of
the price of exotic crosscurrency interest rate derivatives with
pathdependent features via a Partial Differential Equation (PDE)
approach. In particular, we focus on the parallel pricing on Graphics
Processing Unit (GPU) clusters of longdated foreign exchange (FX)
interest rate derivatives, namely PowerReverse DualCurrency (PRDC)
swaps with FX Target Redemption (FXTARN) features under a threefactor
model. Challenges in pricing these derivatives via a PDE approach
arise from the highdimensionality of the model PDE, as well as from
pathdependency of the FXTARN feature. The PDE pricing framework for
FXTARN PRDC swaps 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. Finite difference methods on nonuniform 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 subproblems 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 GPUBased Parallel Algorithm for
Pricing MultiAsset American Options

Authors:  Duy Minh Dang, Christina C. Christara and Kenneth R. Jackson

Date:  March 2011

Pages:  14

Keywords:  American Option, MultiAsset, Penalty Method,
Alternating Direction Implicit Approximate Factorization
(ADIAF), 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 highperformance computational finance, Vol. 24, No. 8, June 2012,
pp. 849866.

File: 
Link to the paper on the SSRN

Abstract:
We develop highlyefficient parallel
Partial Differential Equation (PDE) based pricing methods
on Graphics Processing Units (GPUs) for multiasset 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 GPUbased 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 highorder 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 73113.

File:  ccc/americanOption.pdf

Abstract:
We develop spacetime adaptive and highorder 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 CrankNicolson, are used for the time discretization.
The highorder 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
fourthorder. 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,
multiasset options, multicurrency swaps

Note:  Appeared in the
Canadian Applied Mathematics Quarterly (CAMQ),
Vol. 17, No. 4, 2009, pp. 627660,
which was published in 2011, for a belated 30year 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) timediscretization methods for
solving timedependent 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 GPUbased 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 CrossCurrency 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 crosscurrency interest rate
derivatives under the HullWhite model. In particular, we focus on pricing
longdated 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 CrankNicolson 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 MultiAsset 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, MultiAsset, Penalty Method,
Alternating Direction Implicit Approximate Factorization
(ADIAF), 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, 1319 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 multiasset 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 CrankNicolson, are used for the time discretization. The discrete
nonlinear penalized equations at each timestep are solved using a
penalty iteration. A GPUbased 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 secondorder boundary value problems

Authors:  Christina C. Christara and Guohong Liu 
Date:  September 2009

Pages:  8

Keywords:  sixth order convergence, quartic splines,
secondorder BVP, collocation

Note:  Proceedings of the
9th HERCMA Conference on Computer Mathematics and Applications,
Athens University of Economics and Business,
September 2326, 2009, Athens, Greece, pp. 18.

File:  ccc/paperqrtbvp2.pdf

Abstract:
Collocation methods based on quartic splines are presented for
secondorder twopoint 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 nearboundary
conditions are introduced.
Nonoptimal (fourthorder) and optimal (sixthorder)
quarticspline methods are considered.
The theoretical behavior of the collocation methods
is verified by numerical experiments.
The extension of the methods to twodimensional problems
is briefly considered.
Abstract:
New methods for solving
general linear parabolic partial differential equations (PDEs)
in one space dimension are developed.
The methods combine
quadraticspline collocation for the space discretization and
classical finite differences, such as CrankNicolson,
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 fourthorder 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 15, 2008, Kalamata, Greece, pp. 6267.

File:  ccc/paperqrtbvp4.pdf

Abstract:
We consider the numerical solution
of linear fourthorder 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 HighPerformance 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 912, 2007 Météopole, Toulouse, France, pp. 3537.

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 sixthorder BiQuartic 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:34, 2006, pp. 227257.

File:  SCnonuni.pdf,
DOI

Keywords: 
spline collocation; secondorder twopoint boundary value problem;
error bounds; optimal order of convergence; spline interpolation.

Abstract:
We develop optimal quadratic and cubic spline collocation methods
for solving linear secondorder twopoint boundary value problems
on nonuniform partitions.
To develop optimal nonuniform partition methods,
we use a mapping function from uniform to nonuniform partitions
and develop expansions of the error at the nonuniform 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 nonlinear 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:34, 2006, pp. 259277.

File:  adapt.pdf,
DOI

Keywords: 
spline collocation; secondorder twopoint 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 secondorder twopoint 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 nonuniform points,
placed appropriately to minimize a certain norm of the error.
One adaptive grid technique for cubic spline collocation
is mappingfree 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:  GIAinduced secular variations in the Earth's long wavelength
gravity field: Influence of 3D 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 322327

File:  DOI link ,
journal page

Keywords: 
geopotential harmonics, glacial isostatic adjustment, 3D 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 finitevolume formulation for computing the response
of 3D Maxwell viscoelastic Earth models.
The geometry of the viscosity field is constrained from seismictomographic
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 largescale lateral viscosity variations
of two to three orders of magnitude have a modest, 510%,
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 3D viscoelastic Earth models.
Title:  Glacial isostatic adjustment on 3D Earth models:
a finitevolume 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 421444

File:  DOI link ,
journal page

Keywords: 
crystal motions, finite volumes, glacial isostatic adjustment,
3D viscoelastic Earth models, parallel computation,
domain decomposition

Abstract:
We describe and present results from a finitevolume (FV)
parallel computer code for forward modelling the Maxwell viscoelastic response
of a 3D, selfgravitating, 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 loworder,
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 3D 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 (3D crustal motions, longwavelength gravity anomalies).
These calculations, based on either a simple disc load history
or a global Late Pleistocene ice load reconstruction (ICE3G),
are benchmarked against predictions generated using the traditional
normalmode 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 3D 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 largeamplitude (>1 cm yr^{1}) radial velocities
in zones of peak postglacial 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 3D problem. Specifically, we consider
the impact of mantle viscosity heterogeneity on predictions of
presentday 3D crustal motions in North America. In these tests,
the lateral viscosity variation is constructed, with suitable scaling,
from tomographic images of seismic Swave heterogeneity, and
it is characterized by approximately 2 orders of magnitude (peaktopeak)
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 3D 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 (1D) radial reference model and it far exceeds
observational uncertainties currently being obtained from spacegeodetic
surveying. The relative impact of lateral viscosity variations
on predicted radial motions is significantly smaller.
Title:  Optimal Quadratic Spline Collocation on NonUniform 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; secondorder twopoint 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 secondorder twopoint Boundary Value Problems (BVPs)
discretized on uniform partitions.
In this paper, we extend the optimal QSC methods to nonuniform
partitions, and in a companion paper
we do the same for CSC.
To develop optimal nonuniform partition QSC methods,
we use a mapping function from uniform to nonuniform
partitions and develop expansions of the error at the nonuniform 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^{3j})$ globally, and
$O(h^{4j})$ locally on certain points, for $j \ge 0$.
Title:  Optimal Cubic Spline Collocation on NonUniform 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; secondorder twopoint boundary value problem;
error bounds; optimal order of convergence; adaptive grid;
grid size estimator; error estimator; spline interpolation.

Abstract:
In a recent paper,
nonuniform partition
Quadratic Spline Collocation (QSC) methods
of optimal order of convergence were developed for
secondorder twopoint Boundary Value Problems (BVPs).
In this paper, we develop optimal nonuniform
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 nonuniform partitions.
However, the CSC method is turned into a mappingfree method.
The $j$th derivative of the CSC approximation
is $O(h^{4j})$ globally, for $j \ge 0$, and
$O(h^{5j})$ locally on certain points, for $j > 0$.
The nonuniform 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 mappingfree
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. 175186

File: 
sweqsg02.ps.gz

Keywords:  numerical weather prediction; finite element;
semiLagrangian semiimplicit; Rossby stability; staggered grids

Abstract:
Currently in most global meteorological applications,
the spectral transform method or loworder
finite difference/finite element methods are used.
The spectral transform method, which yields highorder 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
twolevel semiLagrangian semiimplicit 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. 187205

File: 
sweqsc02.ps.gz

Keywords:  numerical weather prediction; finite element;
semiLagrangian semiimplicit; 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 onestep OQSC method, the discrete differential operators are
perturbed to eliminate loworder error terms, and a highorder
approximation is computed using the perturbed operators.
In the twostep OQSC method, a secondorder approximation is generated
first, using the standard formulation, and then a highorder
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 semiLagrangian semiimplicit 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 Ctype grid, the methods also
faithfully capture the group velocity of inertiagravity 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 secondorder 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 secondorder 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 twostep 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 onestep 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, 1316 June 1999.

File:  transpose.99b.ps.gz

Abstract:
Data transposition is required in many numerical applications. When
implemented on a distributedmemory computer, data transposition
requires alltoall 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 shallowwater equations. We
compare the performance of these schemes with that of the Direct
Exchange Algorithm and the MPI alltoall communication routine,
MPI_AllToAll. The numerical experiments were performed on a Cray T3E
computer with 512 processors and on an ethernetconnected 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 wellknown Direct and Ring Exchange
schemes or the MPI_AllToAll routine. We also generalize the Mesh and
Cube Exchange algorithms to a ddimensional 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. 781803.

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 onedimensional secondorder
differential equations.
The rate of convergence of the twogrid method
integrated with a damped Richardson relaxation scheme
as smoother is shown to be faster than 1/2,
independently of the stepsize.
The additive multilevel versions of the
algorithms are also analyzed.
The development of quadratic spline collocation
multigrid methods is extended to twodimensional
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 Vcycle 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. 7189

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 semiiterative
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. DCS135, 1990)

Pages:  34

Note:  A slightly revised version of this paper appeared in BIT, 34:1, 1994, pp. 3361

File:  ccc/cs90135.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 nonoptimal 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 3j )$ global error estimates for the $j$th partial
derivative are obtained for a certain class of problems. Moreover,
$O(h sup 4j )$ 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 CivilComp,
Vol. A, pp. 2334

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 semiiterative
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. 142149
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 35, 1991, Fredericton, NB, Canada, pp. 221230

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 messagepassing 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. 369377

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. 108120, 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 SCPCGQSC method
in its sequential mode is very efficient compared to standard
direct band solvers for the QSC equations.
We have implemented the SCPCGQSC 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 812, 1990, Charleston, SC, U.S.A., pp. 550558

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 CGQSC
(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 68, 1989, Monterey, CA, U.S.A., pp. 12671273.

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 nonshared 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, SaintMalo, France, pp. 375384.

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 collocationcapacitance
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. 175182.

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. 935952.

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^{3j})$ global estimate for the $j$th
derivative of the error is derived.
Further, $O(h^{4j})$ 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^{2j})$ 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. 123155.

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. 179,
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{GolubVanLoan, Parlett, Saad, Stewart, Wilkinson}.
Supervised students' theses
Title: 
Highorder numerical methods for parabolic PDEs
with nonsmooth data: a perspective from option pricing
PDE modelling and computation,

Author:  Dawei Wang

Date:  June 2024

Pages:  124

Note:  Ph.D. thesis, Computer Science Department,
University of Toronto

File: 
dwang24phd.pdf

Abstract:
Highorder methods for solving partial differential
equations (PDEs) enjoy, at least asymptotically, better performance
compared to loworder methods,
when the performance is measured by the ratio of error to computational cost.
However, when there is nonsmoothness in the
solution, as seen in many realworld applications,
it generally requires careful development to solve the
PDEs to a highorder accuracy.
In this thesis, we focus on developing highorder
numerical methods
for solving parabolic PDEs with nonsmooth
data from a perspective of option pricing.
We study separately the convergence behaviors
when the solution contains an unknown free
boundary on which it is nonsmooth (the socalled free
boundary problems), and when the initial conditions are nonsmooth.
For free boundary problems, the finite difference
approximations of derivatives of a nonsmooth
function exhibit degenerated orders of accuracy.
We propose a highorder deferred correction algorithm
combined with penalty iterations, and
show that the order of convergence of the solution can be increased
to fourthorder by solving successively corrected finite difference systems,
where the corrections are derived from the previously computed
lower order solutions, and
applied solely to the righthand side of the linear system.
For handling nonsmooth initial conditions,
while still obtaining highorder time stepping,
we propose to discretize the initial conditions with
appropriate highorder smoothing schemes,
and apply BDF4 time marching initialized with two steps
of an explicit thirdorder RungeKutta (RK3) method and
one step of BDF3 (2RK3BDF3BDF4). From the Fourier analysis
of the discrete system, we prove that, for nonsmooth data,
the loworder errors in the highfrequency domain
are exponentially damped away by BDF steps, while
the persisting loworder errors in the lowfrequency domain
can be eliminated by smoothing techniques, e.g., convolutionbased techniques.
In addition, we derive novel smoothing techniques that
cancel out the loworder terms of the quantization errors
in the Fourier domain arising from discretization.
Furthermore, we show how to flexibly apply smoothings
to general nonsmooth (but piecewise smooth) initial condition discretizations,
not only on uniform, but also on nonuniform grids.
Abundant numerical examples are presented to numerically
support the highorder convergence of the proposed algorithms in the thesis.
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: 
ywchen23phd.pdf

Abstract:
Adjusting derivative prices to take into account default risk
has attracted the attention of several researchers and practitioners,
especially after the 20072008 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 doublepenalty 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
CoxIngersollRoss (CIR) process, while the other party
has constant default risk.
The timedependent PDE derived is of BlackScholes 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 penaltylike 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 meanreversion 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:  rwu21msc.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
HamiltonJacobiBellman equations and introduce new penaltylike 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 secondorder 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 penaltylike methods that we have introduced and
demonstrate that our numerical implementations uphold this property.
Title: 
Numerical Methods for Pricing MultiAsset Options

Author:  Yuwei Chen

Date:  January 2018

Pages:  75

Note:  M.Sc. thesis, Computer Science Department,
University of Toronto

File: 
ywchen18msc.pdf

Abstract:
We consider the pricing of twoasset 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 twoasset options, the basic numerical PDE model
is the twodimensional BlackScholes 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 nonuniform 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 secondorder, for appropriately chosen grid
discretizations.
In the American pricing problem, the discrete penalty method
is applied to the linear complementarity problem involving
the twodimensional BlackScholes PDE and additional constraints.
The convergence of the American Spread put option approximation
computed with penalty iteration remains secondorder, with
the number of penalty iterations per timestep remaining small (23).
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 ChunHo Leung

Date:  June 2017

Pages:  119

Note:  Ph.D. Thesis, Depart. of Computer Science,
University of Toronto, 2017.

File: 
nat17phd.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 nonsmooth 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 nonsmoothness that are of financial interest.
We show that the error of the numerical solution under
CrankNicolsonRannacher 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 nonsmoothness on the grid. We discuss how this positioning
affects the quality of the numerical solution, and the possibility of
an optimal placement of the nonsmoothness 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
timedependent 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: 
mufanLi15basc.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 nonlinear constraint on the BlackScholes 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 CrankNicolson time stepping and nonuniform grids,
the methods are compared in numerical experiments.
The criteria include order of convergence, efficiency, and complexity.
Title: 
Modeling MultiFactor 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: 
crosscurrency interest rate derivatives,
Power Reverse Dual Currency (PRDC) swaps,
Bermudan cancelability, knockout, and FX Target Redemption,
Alternating Direction Implicit (ADI) timestepping methods,
multiGPU platforms, multiasset American options.

Abstract:
This thesis develops efficient modeling frameworks via a Partial
Differential Equation (PDE) approach for multifactor financial
derivatives, with emphasis on threefactor models, and studies
highly efficient implementations of the numerical methods on novel
highperformance computer architectures, with particular focus on Graphics
Processing Units (GPUs) and multiGPU platforms/clusters of GPUs. Two
important classes of multifactor financial instruments are considered:
crosscurrency/foreign exchange (FX) interest rate derivatives and
multiasset options.
For crosscurrency 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 onefactor Gaussian
models for the domestic and foreign interest short rates, and a onefactor
skew model for the spot FX rate results in a timedependent 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 GPUbased 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 multiGPU
platforms/clusters of GPUs to solve each of these independent subproblems
on a separate GPU. Numerical results indicate that the proposed GPUbased
parallel numerical methods are highly efficient and provide significant
increase in performance over CPUbased 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 multiasset options under the BlackScholesMerton framework, with
strong emphasis on multiasset 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 GPUbased 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 GPUbased 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/tyee09msc.pdf

Keywords: 
molecular dynamics, sand, powder, particle simulation, ODE solver, RungeKutta,
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 interparticle 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
generalpurpose granular solver.
Title: 
Multigrid and Spline Collocation Methods on Nonuniform 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
onedimensional secondorder differential equations on
nonuniform 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 nonuniform grids. Multigrid
methods are applied to quadratic and cubic spline collocation equations
arising from the discretization of onedimensional secondorder
differential equations. The convergence analysis of a twogrid
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 secondorder 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 onedimensional secondorder
differential equations discretized on nonuniform grids, and
the conditions under which they are convergent and efficient
are studied. The numerical results confirm our analysis.
Title:  Biquartic Spline Collocation Methods for Fourthorder
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/jingrui08phd.ps.gz

Keywords: 
sixth order of convergence, quartic splines, collocation,
fast Fourier Transform (FFT), GMRES, preconditioning

Abstract:
We present optimal biquartic spline collocation methods for
the numerical solution of fourthorder 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 secondorder approximations, the two
biquartic 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 onestep
(extrapolated) method forces the approximations to satisfy a
perturbed problem, while the threestep (deferredcorrection)
method proceeds in three steps and perturbs the right sides of the
linear equations appropriately.
The threestep biquartic 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 onedimensional problem related to the two
opposite boundaries. The linear system arising from the biquartic
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:  dmdang07msc.pdf

Keywords: 
adaptive mesh selection, monitor function, finite differences, CrankNicolson,
European option, American option, penalty method, projected SOR

Abstract:
We develop spacetime adaptive methods for valuing
American options with strong emphasis on American put options.
We examine the application of adaptive techniques
to the BlackScholes partial differential equation
problem associated with an American put option
in the context of nonuniform secondorder 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 nonuniform grids
and redistribution of the spatial points lead to
linear complementarity problems with $\Mcal$matrices.
The Projected Successive Overrelaxation
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 SecondOrder
TwoPoint Boundary Value ODE Problems

Author:  Guohong Liu

Date:  August 2007

Pages:  130

Note:  M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2007.

File:  ghliu07msc.ps.gz

Keywords: 
sixth order of convergence, quartic splines, collocation,
boundary conditions

Abstract:
Collocation methods based on quartic splines are presented for
secondorder twopoint 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 (onestep) and the deferredcorrection
(twostep) 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. Nonoptimal
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/moqk06msc.pdf

Keywords: 
convertible bond, dividend protection, credit risk,
Conversion Ratio Adjustment, Dividend PassThru, finite differences,
CrankNicolson, Projected Successive OverRelaxation (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 PassThru.
We apply finite difference methods to discretize the PDEs associated
with our models, and study the Projected Successive OverRelaxation (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/tchen05msc.pdf

Keywords: 
semiexplicit semiimplicit method, optimal order of convergence,
CrankNicolson, 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 CrankNicolson,
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:  lucy05msc.pdf

Keywords: 
Coupon payments, call, put, convert, discontinuities,
implicit methods, successive overrelaxation, 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 fixedincome 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 OverRelaxation (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 NonRectangular Domains

Author:  Kit Sun Ng

Date:  June 2005

Pages:  187

Note:  Ph.D. Thesis, Computer Science Dept., Univ. of Toronto, 2005.

File:  ngkit05phd.ps.gz

Keywords: 
quadratic splines, cubic splines, optimal convergence, nonuniform grid,
adaptive methods, gridsize estimator, error estimator,
Lshaped domains, Tshaped 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 nonuniform partitions
to solve onedimensional linear secondorder Boundary Value Problems (BVPs)
that have solutions with regions of high variations, such as boundary layers.
Optimal quadratic spline collocation (QSC) is extended to
nonuniform partitions by using a mapping function from uniform
to nonuniform partitions.
Optimal cubic spline collocation (CSC) on nonuniform 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 nonuniform partitions. However, this method is converted into
a mappingfree method and is integrated with adaptive techniques
that include error and grid size estimators.
Spline collocation is then considered for solving twodimensional
linear secondorder BVPs. We extend optimal spline collocation
to solving BVPs on L and Tshaped 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 onedimensional 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
secondorder BVP into a system of three firstorder 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/zengz05msc.ps.gz

Keywords: 
Shallow water equations, semiLagrangian methods, algebraic elimination,
analytic elimination, systems of PDEs

Abstract:
This thesis describes numerical methods for the solution of
advection equations in onedimensional 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 onedimensional
secondorder differential equations.
The convergence analysis of a twogrid 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
onedimensional shallow water equations (SWEs).
The SWEs are discretized in time with a semiLagrangian
semiimplicit 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 highperformance 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/jingrui04msc.ps.gz

Keywords: 
spline collocation; fourthorder twopoint 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
highperformance 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
onedimensional 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
twodimensional biharmonic Dirichlet problem.
Title:  QuarticSpline Collocation Methods for FourthOrder TwoPoint
Boundary Value Problems

Author:  Ying Zhu

Date:  April 2001

Pages:  73

Note:  M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2001.

File: 
ccc/yz01msc.ps.gz

Abstract:
This thesis presents numerical methods for the solution of
general linear fourthorder 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 quarticspline collocation method is shown
to be second order. Two sixthorder quarticspline collocation methods
are developed and analyzed. They are both based on a high order
perturbation of the differential equation and boundary conditions
operators. The onestep method forces the approximation to satisfy
a perturbed problem, while the threestep 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 quadraticspline collocation matrices.
Numerical results verify the orders of convergence predicted by the
analysis.
Title:  HighOrder 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: 
tam01phd.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 twolevel semiLagrangian
semiimplicit (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
fourthorder in space locally at the nodes and midpoints of the spatial
grids, and thirdorder 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 Cgrids. From an eigenvalue analysis of our
simplified version of the shallow water equations, we conclude that,
when discretized on the Cgrid, our algorithms faithfully capture the
group velocity of inertiagravity waves. Our analysis suggests that
neither the A nor Bgrids 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 Bgrids,
some waves are observed to travel in the wrong direction, whereas, with
the Cgrid, 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/ngkit00msc.ps.gz

Abstract:
We consider Quadratic Spline Collocation (QSC) methods
for solving systems of two linear secondorder 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 ShallowWater 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 shallowwater 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 shallowwater equations on distributed memory
systems. Time integration of the numerical methods is based on a
twotimelevel semiimplicit semiLagrangian scheme. To solve the
Helmholtz problem that arises at each timestep, 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 alltoall 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
ethernetconnected 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 shallowwater 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/constas96msc.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.
Twodimensional and threedimensional 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.