This site provides access to the Technical Reports of the Numerical Analysis and Scientific Computing Group of the Department of Computer Science at the University of Toronto. The reports are listed from newest to oldest.
You may also wish to look at the main page of the Numerical Analysis and Scientific Computing Group.
2017
2016
2015
2014
2013
2012
2011
2010
2009
2008
2004
2001
2000
1999
1998
1997
1996
1995
1994
1993
1992
1991
1990
Title: | A Dimension and Variance Reduction Monte-Carlo Method For Option Pricing Under Jump-Diffusion Models |
Authors: | Duy-Minh Dang, Kenneth R. Jackson and Scott Sues |
Date: | September 2017 (An earlier draft of this paper appeared in October 2015.) |
Pages: | 35 |
Keywords: | conditional Monte Carlo, variance reduction, dimension reduction, partial-integro-differential equations, jump diffusions, fast Fourier transform, normal, double-exponential |
Note: | To appear in the Journal of Applied Mathematical Finance. |
File: | Link to the paper on the SSRN |
We develop a highly efficient MC method for computing plain vanilla European option prices and hedging parameters under a very general jump-diffusion option pricing model which includes stochastic variance and multi-factor Gaussian interest short rate(s). The focus of our MC approach is variance reduction via dimension reduction. More specifically, the option price is expressed as an expectation of a unique solution to a conditional Partial Integro-Differential Equation (PIDE), which is then solved using a Fourier transform technique. Important features of our approach are (i) the analytical tractability of the conditional PIDE is fully determined by that of the Black-Scholes-Merton model augmented with the same jump component as in our model, and (ii) the variances associated with all the interest rate factors are completely removed when evaluating the expectation via iterated conditioning applied to only the Brownian motion associated with the variance factor. For certain cases when numerical methods are either needed or preferred, we propose a discrete fast Fourier transform method to numerically solve the conditional PIDE efficiently. Our method can also effectively compute hedging parameters. Numerical results show that the proposed method is highly efficient.
Title: | Correlated Multivariate Poisson Processes and Extreme Measures |
Authors: | Michael Chiu, Kenneth R. Jackson and Alexander Kreinin |
Date: | June 2017 |
Pages: | 13 |
Note: | Submitted to the Journal of Model Assisted Statistics and Applications (MASA). |
File: | Correlated_Poisson_Processes_1.pdf |
Multivariate Poisson processes have many important applications in Insurance, Finance, and many other areas of Applied Probability. In this paper we study the backward simulation approach to modelling multivariate Poisson processes and analyze the connection to the extreme measures describing the joint distribution of the processes at the terminal simulation time.
Title: | Mixing LSMC and PDE Methods to Price Bermudan Option |
Authors: | David Farahany, Kenneth R. Jackson and Sebastian Jaimungal |
Date: | November 2016; revised August 2017 |
Pages: | 33 |
Note: | Work In Progress |
File: | Link to the paper on the SSRN |
We develop a mixed least square Monte Carlo-partial differential equation (LSMC-PDE) method for pricing Bermudan style options with assets whose dynamics are driven by latent variables. The algorithm is formulated for an arbitrary number of assets and volatility processes, and its probabilistic convergence is established. Our numerical examples focus on the Heston model and we compare our hybrid algorithm with a full LSMC approach. Using Fourier methods we are able to derive efficient FFT based solutions, and we demonstrate that our algorithm greatly reduces the variance in the computed prices and optimal exercise boundaries. We also compare the early exercise boundaries and prices computed by our hybrid algorithm with those produced by finite difference methods and find excellent agreement.
Title: | Numerical Methods for Computing Sensitivities for ODEs and DDEs |
Author: | J. Calver and W.H. Enright |
Date: | September 2016 |
Pages: | 17 |
Keywords: | Ordinary differential equations, delay differential equations, adjoint method, variational equations, sensitivities, parameter estimation |
Note: | A revised version of this paper will appear in Numerical Algorithms |
File: | adjoint_paperNSERC.pdf |
We investigate the performance of the adjoint approach and the variational approach for computing the sensitivities of the least squares ob- jective function commonly used when fitting models to observations. We note that the discrete nature of the objective function makes the cost of the adjoint approach for computing the sensitivities dependent on the number of observa- tions. In the case of ODEs, this dependence is due to having to interrupt the computation at each observation point during numerical solution of the adjoint equations. Each observation introduces a jump discontinuity in the solution of the adjoint differential equations. These discontinuities are propagated in the case of DDEs, making the performance of the adjoint approach even more sensitive to the number of observations for DDEs. We quantify this cost and suggest ways to make the adjoint approach scale better with the number of observations. In numerical experiments, we compare the adjoint approach with the variational approach for computing the sensitivities.
Title: | A Neural Network Approach to Efficient Valuation of Large VA Portfolios |
Author: | Seyed Amir Hejazi |
Date: | August 2016 |
Pages: | 162 |
Keywords: | Variable annuity, Spatial interpolation, Neural network, Portfolio valuation, Solvency Capital Requirement (SCR) |
Note: | PhD Thesis |
File: | Seyed_Amir_Hejazi_2016_PhD_Thesis.pdf |
Variable annuity (VA) products expose insurance companies to considerable risk because of the guarantees they provide to buyers of these products. Managing and hedging the risks associated with VA products requires intraday valuation of key risk metrics for these products. The complex structure of VA products and computational complexity of their accurate evaluation has compelled insurance companies to adopt Monte Carlo (MC) simulations to value their large portfolios of VA products. Because the MC simulations are computationally demanding, especially for intraday valuations, insurance companies need more efficient valuation techniques.
Existing academic methodologies focus on fair valuation of a single VA contract, exploiting ideas in option theory and regression. In most cases, the computational complexity of these methods surpasses the computational requirements of MC simulations. Recently, a framework based on Kriging has been proposed that can significantly decrease the computational complexity of MC simulation. Kriging methods are an important class of spatial interpolation techniques. In this thesis, we study the performance of prominent traditional spatial interpolation techniques. Our study shows that traditional interpolation techniques require the definition of a distance function that can significantly impact their accuracy. Moreover, none of the traditional spatial interpolation techniques provide all of the key properties of accuracy, efficiency, and granularity. Therefore, in this thesis, we present a neural network approach for the spatial interpolation framework that affords an efficient way to find an effective distance function. The proposed approach is accurate, efficient, and provides an accurate granular view of the input portfolio. Our numerical experiments illustrate the superiority of the performance of the proposed neural network approach in estimation of the delta value and also the solvency capital requirement for large portfolios of VA products compared to the traditional spatial interpolation schemes and MC simulations.
Title: | Limiting Deviation Method for Coupling Codes with Adaptive Window-Size Control using Digital Filters |
Author: | Rohan Palaniappan |
Date: | July 2015 |
Pages: | 75 |
Note: | MSc Research Paper, Computer Science Department, University of Toronto |
File: | rohan-palaniappan-msc-research-paper.pdf |
Multi-physics systems model multiple simultaneous physical phenomena that affect one another. Simulating such systems can be done using the Operator Splitting (OS) method, wherein separate simulations are run concurrently for each phenomenon and data is exchanged between the simulations periodically. In this work, we explore, in the context of the OS method, what we call the Limiting Deviation (LD) and Limiting Deviation with Interior Check (LDIC) methods, which indirectly limit the global solution error in an unconventional way. We test the LD and LDIC methods using "canned data" for five test cases using thirteen different adaptive step-size controllers. We identify the best controllers for two experiments, point out undesired behaviour in certain controllers and examine the effect of two key parameters. Finally, we discuss implementation options. While our results are quite positive, we consider them preliminary, because the tests do not factor-in the coupling between the separate simulations.
Title: | Efficient and Accurate Numerical PDE Methods for Pricing Financial Derivatives, |
Author: | Mufan (Bill) Li |
Date: | April 2015 |
Pages: | 35 |
Note: | B.A.Sc. Thesis, Depart. of Electrical and Computer Eng., Univ. of Toronto, 2015. |
File: | mufanLi-15-basc.pdf |
Keywords: | American options, (implicit) penalty iteration, operator splitting (explicity penalty) method, convergence, accuracy, complexity, efficiency, Greek calculation. |
Title: | Efficient Valuation of SCR via a Neural Network Approach |
Author: | Seyed Amir Hejazi and Kenneth R. Jackson |
Date: | First published September 2015; revised October 2016. |
Pages: | 29 |
Keywords: | Variable annuity, Spatial interpolation, Neural network, Portfolio valuation, Solvency Capital Requirement (SCR) |
Note: | Accepted for publication in
Journal of Computational and Applied Mathematics. An earlier version of this paper can be found here. |
File: | CAM-Paper3.pdf |
As part of the new regulatory framework of Solvency II, introduced by the European Union, insurance companies are required to monitor their solvency by computing a key risk metric called the Solvency Capital Requirement (SCR). The official description of the SCR is not rigorous and has lead researchers to develop their own mathematical frameworks for calculation of the SCR. These frameworks are complex and are difficult to implement. Recently, Bauer et al. suggested a nested Monte Carlo (MC) simulation framework to calculate the SCR. But the proposed MC framework is computationally expensive even for a simple insurance product. In this paper, we propose incorporating a neural network approach into the nested simulation framework to significantly reduce the computational complexity in the calculation. We study the performance of our neural network approach in estimating the SCR for a large portfolio of an important class of insurance products called Variable Annuities (VAs). Our experiments show that the proposed neural network approach is both efficient and accurate.
Title: | A Neural Network Approch to Efficient Valuation of Large Portfolios of Variable Annuities |
Author: | Seyed Amir Hejazi and Kenneth R. Jackson |
Date: | May 2016 (earlier drafts September 2015 and January 2016) |
Pages: | 37 |
Keywords: | Variable annuity, Spatial interpolation, Neural Network, Portfolio valuation |
Note: | To appear in Insurance: Mathematics and Economics |
File: | IME-Paper2b.pdf
(Earlier versions of this paper can be found in draft1 and draft2.) |
Managing and hedging the risks associated with Variable Annuity (VA) products require intraday valuation of key risk metrics for these products. The complex structure of VA products and computational complexity of their accurate evaluation have compelled insurance companies to adopt Monte Carlo (MC) simulations to value their large portfolios of VA products. Because the MC simulations are computationally demanding, especially for intraday valuations, insurance companies need more efficient valuation techniques. Recently, a framework based on traditional spatial interpolation techniques has been proposed that can significantly decrease the computational complexity of MC simulation (Gan and Lin, 2015). However, traditional interpolation techniques require the definition of a distance function that can significantly impact their accuracy. Moreover, none of the traditional spatial interpolation techniques provide all of the key properties of accuracy, efficiency, and granularity (Hejazi et al., 2015). In this paper, we present a neural network approach for the spatial interpolation framework that affords an efficient way to find an effective distance function. The proposed approach is accurate, efficient, and provides an accurate granular view of the input portfolio. Our numerical experiments illustrate the superiority of the performance of the proposed neural network approach compared to the traditional spatial interpolation schemes.
Title: | Computational Modeling and Analysis of Complex Muscle Architecture |
Author: | Dongwoon Lee |
Date: | August 2015 |
Pages: | 94 |
Keywords: | skeletal muscle; muscle architecture; cadaveric specimen; PCSA; pennation angle |
Note: | PhD Thesis |
File: | Lee_Dongwoon_PhDThesis.pdf |
Muscle architecture is a primary determinant of the muscle function associated with body movement. An assessment of muscle architecture is therefore of great importance, not only for investigating anatomical aspects of muscle but also for predicting its functional capacity. Most muscles have a variable complexity in their architectures, making it challenging to accurately assess them. Previous cadaveric approaches only take into account limited portions of architecture. On the other hand, conventional radiological approaches, such as ultrasonography and MRI, examine two-dimensional projected images. Neither of these approaches provides a thorough understanding of the entire muscle architecture. This may lead to under- or over-estimation of architectural parameters that are significant for both clinical and computational studies. Therefore, the purpose of this thesis is to develop a computational modeling approach to facilitate quantification and reconstruction of complex muscle architecture. Cadaveric specimen data are used to investigate muscle architecture and to reconstruct accurate models. Associated geometric complexity and variation are carefully examined to yield consistent estimation of architectural parameters. This method demonstrates robustness against non-uniformity in the data and consistency over various types of muscle architecture; less than 10% error in PCSA estimation. By incorporating ultrasonographic assessment, this method is extended to approximate muscle architecture in living subjects, which enables estimation of PCSA for in vivo muscle in a more consistent manner. Validation experiments demonstrate 0.4% - 8.4% estimation errors between the original architecture and its approximation, depending on the anatomical complexity, which provides a practical insight into the quantification of PCSA for in vivo muscle.
Title: | New Approaches to Importance Sampling for Portfolio Credit Risk Valuation |
Author: | Zhe (Robert) Wang |
Date: | July 2015 |
Pages: | 57 |
Keywords: | portfolio credit risk; Gaussian copula; Monte Carlo simulation; importance sampling; variance reduction |
Note: | MSc Research Paper |
File: | Zhe.Wang.MSC.Thesis.pdf |
Portfolio credit risk based on the Gaussian Copula model has been widely studied and generally evaluated through Monte Carlo simulations. The two-level structure, namely systematic factors and individual factors, complicates the problem in a way that traditional variance reduction techniques become very hard to apply. Glasserman and Li proposed a two-level importance sampling approach to tackle a simplified binary credit states problem. The inner level was approximated through a conditional importance sampling approach using an exponential twisting technique. In this research project, we propose an alternative importance sampling approach which uses the Central Limit Theorem for the inner level. Our approach can be easily generalized to multi-credit states. Based on this approximation, we then propose two novel approaches motivated from research in machine learning. Instead of finding the importance sampler through an optimization problem, the first approach approximates the zero variance function by learning from the samples which are generated from Markov Chain Monte Carlo. The second approach treats the problem as a Bayesian inference problem and evaluates the tail probability through Bayesian Monte Carlo. Compared to Glasserman and Li's method, numerical results show that these two new approaches have advantages in both accuracy and speed and are also more easily adapted to other applications.
Title: | A Practical Approach to in vivo Quantification of Physiological Cross-Sectional Area for Human Skeletal Muscle |
Authors: | Dongwoon Lee, Soo Kim, Ken Jackson, Eugene Fiume and Anne Agur |
Date: | May 2015 |
Pages: | 9 |
Keywords: | skeletal muscle; muscle architecture; PCSA; ultrasonography; digitization |
Note: | Submitted to the Journal of Biomechanics |
File: | Lee.JBM.2015.pdf |
Physiological cross-sectional area (PCSA) is an important property used to predict the maximal force capacity of skeletal muscle. A common approach to estimate PCSA uses an algebraic formula based on muscle volume (MV), fascicle length (FL) and pennation angle (PA) that are measured by magnetic resonance imaging (MRI) and ultrasonographic assessments. Due to the limited capability of assessing architecturally complex muscles with these imaging modalities, the accuracy of measurements and ultimately PCSA estimation is compromised. On the other hand, cadaveric modeling provides a more accurate quantification by effectively dealing with the variable complexity of the muscle but it may not be straightforward to directly apply to in vivo studies. Thus, the purpose of our study is to provide a practical approach to PCSA estimation for in vivo muscle by integrating both cadaveric and ultrasound data. The muscle architecture in vivo is approximated by fitting 3D cadaveric data onto 2D ultrasound data of living subjects. The fitted architectural data is used for PCSA quantification. Validation experiments based on synthetic muscle and cadaveric data, respectively, demonstrate 0.4 — 8.4% errors between original architecture and its approximation, depending on the anatomical complexity. Furthermore, it is shown that, despite the large inter-subject variability of cadaveric data (standard deviation: ±153.2 mm^{2}), their transformation toward 2D ultrasound data consistently yields a narrow distribution of PCSA estimation (standard deviation: ±24.6 ~ ±35.7 mm^{2}), which provides a practical insight into accurate quantification of PCSA for in vivo muscle.
Title: | A Spatial Interpolation Framework for Efficient Valuation of Large Portfolios of Variable Annuities |
Authors: | Seyed Amir Hejazi, Kenneth R. Jackson and Guojun Gan |
Date: | April 2017 |
Pages: | 27 |
Keywords: | Variable annuity, Spatial interpolation, Kriging, Inverse distance weighting, Radial basis function, Portfolio valuation |
Note: | To appear in the special issue on Computational Finance and Insurance of the Journal of Quantitative Finance and Economics. |
File: | IME-Paper1.pdf |
Variable Annuity (VA) products expose insurance companies to considerable risk because of the guarantees they provide to buyers of these products. Managing and hedging these risks requires insurers to find the value of key risk metrics for a large portfolio of VA products. In practice, many companies rely on nested Monte Carlo (MC) simulations to find key risk metrics. MC simulations are computationally demanding, forcing insurance companies to invest hundreds of thousands of dollars in computational infrastructure per year. Moreover, existing academic methodologies are focused on fair valuation of a single VA contract, exploiting ideas in option theory and regression. In most cases, the computational complexity of these methods surpasses the computational requirements of MC simulations. Therefore, academic methodologies cannot scale well to large portfolios of VA contracts. In this paper, we present a framework for valuing such portfolios based on spatial interpolation. We provide a comprehensive study of this framework with existing interpolation schemes. Our numerical results show superior performance, in terms of both computational efficiency and accuracy, for these methods compared to nested MC simulations. We also present insights into the challenge of finding an effective interpolation scheme in this framework, and suggest guidelines that help us build a fully automated scheme that is efficient and accurate.
Title: | Dimension and Variance Reduction for Monte Carlo Methods for High-Dimensional Models in Finance |
Authors: | Duy-Minh Dang, Kenneth R. Jackson and Mohammadreza Mohammadi |
Date: | January 2015 |
Pages: | 21 |
Keywords: | conditional Monte Carlo, variance reduction, dimension reduction, cross-currency, Fourier transform, partial differential equations |
JEL Classification: | E40, E43, G12, G13, C61, C63 |
File: | Link to the paper on the SSRN |
One-way coupling often occurs in multi-dimensional models in finance. In this paper, we present a dimension reduction technique for Monte Carlo (MC) methods, referred to as drMC, that exploits this structure for pricing plain-vanilla European options under a N-dimensional one-way coupled model, where N is arbitrary. The dimension reduction also often produces a significant variance reduction.
The drMC method is a dimension reduction technique built upon (i) the conditional MC technique applied to one dimension and (ii) the derivation of a closed-form solution for the conditional Partial Differential Equation (PDE) that arises via Fourier transforms. In the drMC approach, the option price can be computed simply by taking the expectation of this closed-form solution. Hence, the approach results in a powerful dimension reduction from N to one, which often results in a significant variance reduction as well, since the variance associated with the other (N-1) factors in the original model are completely removed from the drMC simulation. Moreover, under the drMC framework, hedging parameters, or Greeks, can be computed in a much more efficient way than in traditional MC techniques.
A variance reduction analysis of the method is presented and numerical results illustrating the method's efficiency are provided.
Title: | Improving a Cross Entropy Approach to Parameter Estimation for ODEs and DDEs |
Author: | Jonathan Calver |
Date: | February 2014 |
Pages: | 30 |
Note: | MSc Thesis |
File: | j_calver_MSc_paper.pdf |
We investigate and extend the cross entropy (CE) approach for parameter estimation for ODE and DDE models introduced in [27]. Software is developed to allow models to be easily speciﬁed for use with an existing package for working with DDEs. Our software implements CE and is used to explore potential improvements to the method. A 3-stage optimization procedure is presented. First, an inexpensive initial guess for some of the parameters is obtained to reduce the size of the search space. Second, a global search phase is performed to obtain an estimate for the parameters that is ‘close’ to the optimal values. Third, a local optimizer, that can converge super-linearly, is used to obtain the ﬁnal estimates for the parameters. The performance of this 3-stage procedure is demonstrated by estimating parameters on several test problems from the literature. Additionally, we look at ways to reduce the computational cost of simulating the ODEs and DDEs during the estimation procedure.
Title: | A Survey of Numerical Methods for Lévy Markets |
Author: | Michael Chiu |
Date: | August 2014 |
Pages: | 29 |
Note: | M. Eng. Project, MIE Dept., University of Toronto |
File: | Michael.Chiu.MEng.Project.pdf |
The modeling of financial markets by Lévy proceses has become an active area of research during recent years. This has motivated an equal amount of, if not more, research activity into the necessary numerical methods. Due to the large body of work in this area, we focus our survey on fast numerical methods for Lévy markets. Particular emphasis is placed on grid-based methods.
Title: | A Three-Dimensional Approach to Pennation Angle Estimation for Human Skeletal Muscle |
Authors: | Dongwoon Lee, Zhi Li, Qazi Zain Sohail, Ken Jackson, Eugene Fiume and Anne Agur |
Date: | March 2014 (earlier draft July 2013) |
Pages: | 15 |
Keywords: | pennation angle; line of action; ultrasonography; digitization; cadaveric specimen |
Note: | Published in the Journal of
Computer Methods in Biomechanics and Biomedical Engineering
An earlier version of this paper can be found here. |
File: | Lee.PA.2014.pdf |
Pennation angle (PA) is an important property of human skeletal muscle that plays a significant role in determining the force contribution of fascicles to skeletal movement. Two-dimensional (2D) ultrasonography is the most common approach to measure PA. However, in principle, it is challenging to infer knowledge of three-dimensional (3D) architecture from 2D assessment. Furthermore, architectural complexity and variation impose more difficulties on reliable and consistent quantification of PA. Thus, the purpose of our study is to provide accurate insight into the correspondence between 2D assessment and the underlying 3D architecture. To this end, a 3D method was developed to directly quantify PA based on 3D architectural data that were acquired from cadaveric specimens through dissection and digitization. Those data were then assessed two-dimensionally by simulating ultrasound imaging. To achieve consistency over intermuscular variation, our proposed 3D method is based on the geometric analysis of fascicle attachment. Comparative results show a wide range of differences (1.1% - 47.1%) between 2D and 3D measurements. That is, ultrasound can under- or over-estimate PA, depending on architecture.
Title: | A Highly Efficient Implementation on GPU Clusters of PDE-Based Pricing Methods for Path-Dependent Foreign Exchange Interest Rate Derivatives |
Authors: | Duy Minh Dang, Christina Christara and Kenneth R. Jackson |
Date: | Submitted March 2013, revised April 2013. |
Pages: | 17 |
Keywords: | Power Reverse Dual Currency (PRDC) Swaps, Bermudan Cancelable, Partial Differential Equation (PDE), Alternating Direction Implicit (ADI), Finite Differences, Graphics Processing Units (GPUs), GPU Clusters, MPI, Parallel Computing |
Note: | Published in Springer Lecture Notes in Computer Science, Volume 7975, 2013, pp.~107--126, Proceedings of The 13th International Conference on Computational Science and Its Applications (ICCSA 2013), Ho Chi Minh City, Vietnam, 24--27 June 2013. |
File: | Link to the paper on the SSRN |
We present a highly efficient parallelization of the computation of the price of exotic cross-currency interest rate derivatives with path-dependent features via a Partial Differential Equation (PDE) approach. In particular, we focus on the parallel pricing on Graphics Processing Unit (GPU) clusters of long-dated foreign exchange (FX) interest rate derivatives, namely Power-Reverse Dual-Currency (PRDC) swaps with FX Target Redemption (FX-TARN) features under a three-factor model. Challenges in pricing these derivatives via a PDE approach arise from the high-dimensionality of the model PDE, as well as from path-dependency of the FX-TARN feature. The PDE pricing framework for FX-TARN PRDC swaps is based on partitioning the pricing problem into several independent pricing sub-problems over each time period of the swap's tenor structure, with possible communication at the end of the time period. Finite difference methods on non-uniform grids are used for the spatial discretization of the PDE, and the Alternating Direction Implicit (ADI) technique is employed for the time discretization. Our implementation of the pricing procedure on a GPU cluster involves (i) efficiently solving each independent sub-problems on a GPU via a parallelization of the ADI timestepping technique, and (ii) utilizing MPI for the communication between pricing processes at the end of the time period of the swap's tenor structure. Numerical results showing the efficiency of the parallel methods are provided.
Title: | New Greedy Heuristics for Approximating Set Cover and Set Packing |
Author: | David Kordalewski |
Date: | April 2013 |
Pages: | 58 |
Note: | MSc Thesis, Computer Science Dept., Univ. of Toronto, 2013. |
File: | Kordalewski.MSc.Thesis.pdf |
The Set Cover problem (SCP) and Set Packing problem (SPP) are standard NP-hard combinatorial optimization problems. Their decision problem versions are shown to be NP-Complete in Karp's 1972 paper. We specify a rough guide to constructing approximation heuristics that may have widespread applications and apply it to devise greedy approximation algorithms for SCP and SPP, where the selection heuristic is a variation of that in the standard greedy approximation algorithm. Our technique involves assigning to each input set a valuation and then selecting, in each round, the set whose valuation is highest. We prove that the technique we use for determining a valuation of the input sets yields a unique value for all Set Cover instances. For both SCP and SPP we give experimental evidence that the valuations we specify are unique and can be computed to high precision quickly by an iterative algorithm. Others have experimented with testing the observed approximation ratio of various algorithms over a variety of randomly generated instances, and we have extensive experimental evidence to show the quality of the new algorithm relative to greedy heuristics in common use. Our algorithms are somewhat more computationally intensive than the standard heuristics, though they are still practical for large instances. We discuss some ways to speed up our algorithms that do not significantly distort their effectiveness in practice on random instances.
Title: | Iterative Reconstruction Algorithms for Polyenergetic X-Ray Computerized Tomography |
Author: | Nargol Rezvani |
Date: | August 2012 |
Pages: | 127 |
Note: | PhD Thesis, Computer Science Dept., Univ. of Toronto, 2012. |
File: | Nargol.Rezvani.PhD.Thesis.pdf |
A reconstruction algorithm in computerized tomography is a procedure for reconstructing the attenuation coefficientscient, a real-valued function associated with the object of interest, from the measured projection data. Generally speaking, reconstruction algorithms in CT fall into two categories: direct, e.g., filtered back-projection (FBP), or iterative. In this thesis, we discuss a new fast matrix-free iterative reconstruction method based on a polyenergetic model.
While most modern x-ray CT scanners rely on the well-known filtered back-projection algorithm, the corresponding reconstructions can be corrupted by beam hardening artifacts. These artifacts arise from the unrealistic physical assumption of monoenergetic x-ray beams. In this thesis, to compensate, we use an alternative model that accounts for differential absorption of polyenergetic x-ray photons and discretize it directly. We do not assume any prior knowledge about the physical properties of the scanned object. We study and implement different solvers and nonlinear unconstrained optimization methods, such as a Newton-like method and an extension of the Levenberg-Marquardt-Fletcher algorithm. We explain how we can use the structure of the Radon matrix and the properties of FBP to make our method matrix-free and fast. Finally, we discuss how we regularize our problem by applying different regularization methods, such as Tikhonov and regularization in the 1-norm. We present numerical reconstructions based on the associated nonlinear discrete formulation incorporating various iterative optimization methods.
Title: | An Efficient Numerical PDE Approach for Pricing Foreign Exchange Interest Rate Hybrid Derivatives |
Authors: | Duy Minh Dang, Christina Christara, Kenneth R. Jackson and Asif Lakhany |
Date: | March 2012; revised May 2013 |
Pages: | 36 |
Keywords: | Power-Reverse Dual-Currency (PRDC) swaps, Target Redemption (TARN), knockout, Partial Differential Equation (PDE), finite differences, non-uniform grids. |
Note: | Accepted for publication by the Journal of Computational Finance, Vol 18, Issue 4, 2015, pp 39-93. |
File: | Link to the paper on the SSRN |
We discuss efficient pricing methods via a Partial Differential Equation (PDE) approach for long dated foreign exchange (FX) interest rate hybrids under a three-factor multi-currency pricing model with FX volatility skew. The emphasis of the paper is on Power-Reverse Dual-Currency (PRDC) swaps with popular exotic features, namely knockout and FX Target Redemption (FX-TARN). Challenges in pricing these derivatives via a PDE approach arise from the high-dimensionality of the model PDE, as well as from the complexities in handling the exotic features, especially in the case of the FX-TARN provision, due to its path-dependency. Our proposed PDE pricing framework for FX-TARN PRDC swaps is based on partitioning the pricing problem into several independent pricing sub-problems over each time period of the swap's tenor structure, with possible communication at the end of the time period. Each of these pricing sub-problems can be viewed as equivalent to a knockout PRDC swap with a known time-dependent barrier, and requires a solution of the model PDE, which, in our case, is a time-dependent parabolic PDE in three space dimensions. Finite difference schemes on non-uniform grids are used for the spatial discretization of the model PDE, and the Alternating Direction Implicit (ADI) timestepping methods are employed for its time discretization. Numerical examples illustrating the convergence properties and efficiency of the numerical methods are provided.
Title: | Accurate First-Order Sensitivity Analysis for Delay Differential Equations: Part1: The Forward Approach |
Authors: | H. ZivariPiran and W.H. Enright |
Date: | 2011 |
Pages: | 22 |
Note: | A revised version of this paper will appear in SIAM J. Sci. Comp., 2012. |
File: | ZE_SISC.pdf |
In this paper, we derive an equation governing the dynamics of first-order forward sensitivities for a general system of parametric neutral delay differential equations (NDDEs). We also derive a formula which identifies the size of jumps that appear at discontinuity points when the sensitivity equations are integrated. The formula leads to an algorithm which can compute sensitivities for various types of parameters very accurately and efficiently.
Title: | Improved Implementation of Multiple Shooting for BVPs |
Authors: | Weidong Zhang |
Date: | January, 2012 |
Pages: | 26 |
Note: | MSc. Thesis, Computer Science Dept., Univ. of Toronto, 2012. |
File: | Weidong_Zhang_Thesis.pdf |
Boundary value problems arise in many applications, and shooting methods are one approach to approximate the solution of such problems.A Shooting method transforms a boundary value problem into a sequence of initial value problems,and takes the advantage of the speed and adaptivity of initial value problem solvers. The implementation of continuous Runge-Kutta methods with defect control for initial value problems gives efficient and reliable solutions. In this paper, we design and implement a boundary value solver that is based on a shooting method using a continuous Runge-Kutta method to solve the associated initial value problems.Numerical tests on a selection of problems show that this approach achieves better performance than another widely used existing shooting method.
Title: | Reducing the Uncertainty when Approximating the Solution of ODEs |
Authors: | W.H. Enright |
Date: | 2012 |
Pages: | 13 |
Note: | A revised version of this paper will appear in proceedings of IFIP WG 2.5 Conference on Uncertainty Quantification, Boulder, August 2011. |
File: | E_Boulder.pdf |
One can reduce the uncertainty in the quality of an approximate solution of an ordinary differential equation (ODE) by implementing methods which have a rigorous error control strategy and which deliver an approximate solution that is much more likely to satisfy the expectations of the user. We have developed such a class of ODE methods as well as a collection of software tools that will deliver a piecewise polynomial as the approximate solution and facilitate the investigation of various aspects of the problem that are often of as much interest as the approximate solution itself. We will introduce measures that can be used to quantify the reliability of an approximate solution and discuss how one can implement methods that, at some extra cost, can produce very reliable approximate solutions and therefore significantly reduce the uncertainty in the computed results.
Title: | Superconvergent Interpolants for Collocation Methods applied to Volterra Integro-differential Equations with Delay |
Authors: | M. Shakourifar and W.H. Enright |
Date: | 2012 |
Pages: | 16 |
Note: | A revised version of this paper will appear in BIT 2012. |
File: | SE_BIT.pdf |
Standard software based on the collocation method for differential equations delivers a continuous approximation (called the collocation solution) which augments the high order discrete approximate solution that is provided at mesh points. This continuous approximation is less accurate than the discrete approximation. For 'non-standard' Volterra integro-differential equations with constant delay, that often arise in modeling predator-prey systems in Ecology, the collocation solution is $C^0$ continuous. The accuracy is $O(h^{s+1})$ at off-mesh points and $O(h^{2 s})$ at mesh points where $s$ is the number of Gauss points used per subinterval and $h$ refers to the stepsize. We will show how to construct $C^1$ interpolants with an accuracy at off-mesh points and mesh points of the same order ($2 s$). This implies that even for coarse mesh selections we achieve an accurate and smooth approximate solution. Specific schemes are presented for $s = 2, 3$, and numerical results demonstrate the effectiveness of the new interpolants.
Title: | Robust Estimation of PCSA and Geometric Reconstruction for Human Skeletal Muscle |
Authors: | Dongwoon Lee, Kajeandra Ravichandiran, Ken Jackson, Eugene Fiume and Anne Agur |
Date: | January, 2012 |
Pages: | 10 |
Keywords: | Skeletal muscle; Muscle Architecture; Computational Geometry; Digitization |
Note: | To appear in the Journal of Biomechanics under a slightly different title: "Robust estimation of physiological cross-sectional area and geometric reconstruction for human skeletal muscle". |
File: | Lee.JBM.2012.pdf |
Understanding muscle architecture is crucial to determining the mechanical function of muscle during body movements, because architectural parameters directly correspond to muscle performance. Accurate parameters are thus essential for reliable simulation. Human cadaveric muscle specimen data provides the anatomical detail needed for in-depth understanding of muscle and accurate parameter estimation. However, as muscle generally has non-uniform architecture, parameter estimation, specifically, physiological cross-sectional area (PCSA), is rarely straightforward. To deal ectively with this non-uniformity, we propose a geometric approach in which a polygon is sought to best approximate the cross-sectional area of each fascicle by accounting for its three-dimensional trajectory and arrangement in the muscle. Those polygons are then aggregated to determine PCSA and volume of muscle. Experiments are run using both synthetic data and muscle specimen data. From comparison of PCSA using synthetic data, we conclude that the proposed method enhances the robustness of PCSA estimation against variation in muscle architecture. Furthermore, we suggest reconstruction methods to extract 3D muscle geometry directly from fascicle data and estimated parameters using the level set method.
Title: | Parameter Estimation for ODEs using a Cross-Entropy Approach |
Authors: | Bo Wang |
Date: | January, 2012 |
Pages: | 34 |
Note: | MSc. Thesis, Computer Science Dept., Univ. of Toronto, 2012. |
File: | Bo_Wang_Thesis.pdf |
Parameter Estimation for ODEs and DDEs is an important topic in numerical analysis. In this paper, we present a novel approach to address this inverse problem. Cross-entropy algorithms are general a lgorithm which can be applied to solve global optimization problems. The main steps of cross-entrop y methods are first to generate a set of trial samples from a certain distribution, then to update the distribution based on these generated sample trials. To overcome the prohibitive computation o f standard cross-entropy algorithms, we develop a modification combining local search techniques. T he modified cross-entropy algorithm can speed the convergence rate and improve the accuracy simulta neously. Two different coding schemes (continuous coding and discrete coding) are also introduced. Continuou s coding uses a truncated multivariate Gaussian to generate trial samples, while discrete coding re duces the search space to a finite (but dense) subset of the feasible parameter values and uses a B ernoulli distribution to generate the trial samples (which are fixed point approximation of the act ual parameters) . Extensive numerical and real experiments are conducted to illustrate the power an d advantages of the proposed methods. Compared to other existing state-of-the-art approaches on som e benchmark problems for parameter estimation, our methods have three main advantages: 1) They are robust to noise in the data to be fitted; 2) They are not sensitive to the number of observation po ints (in contrast to most existing approaches) ; 3) The modified versions exhibit faster convergenc e than the original versions, thus they are more efficient, without sacrificing accuracy.
Title: | Modeling Multi-Factor Financial Derivatives by a Partial Differential Equation Approach with Efficient Implementation on Graphics Processing Units |
Author: | Duy Minh Dang |
Date: | August 2011 |
Pages: | 229 |
Note: | PhD Thesis, Computer Science Dept., Univ. of Toronto, 2011. |
File: | Duy_Minh_Dang_PhD_Thesis.pdf |
This thesis develops efficient modeling frameworks via a Partial Differential Equation (PDE) approach for multi-factor financial derivatives, with emphasis on three-factor models, and studies highly efficient implementations of the numerical methods on novel high-performance computer architectures, with particular focus on Graphics Processing Units (GPUs) and multi-GPU platforms/clusters of GPUs. Two important classes of multi-factor financial instruments are considered: cross-currency/foreign exchange (FX) interest rate derivatives and multi-asset options.
For cross-currency interest rate derivatives, the focus of the thesis is on Power Reverse Dual Currency (PRDC) swaps with three of the most popular exotic features, namely Bermudan cancelability, knockout, and FX Target Redemption. The modeling of PRDC swaps using one-factor Gaussian models for the domestic and foreign interest short rates, and a one-factor skew model for the spot FX rate results in a time-dependent parabolic PDE in three space dimensions. Our proposed PDE pricing framework is based on partitioning the pricing problem into several independent pricing subproblems over each time period of the swap's tenor structure, with possible communication at the end of the time period. Each of these subproblems requires a solution of the model PDE. We then develop a highly efficient GPU-based parallelization of the Alternating Direction Implicit (ADI) timestepping methods for solving the model PDE. To further handle the substantially increased computational requirements due to the exotic features, we extend the pricing procedures to multi-GPU platforms/clusters of GPUs to solve each of these independent subproblems on a separate GPU. Numerical results indicate that the proposed GPU-based parallel numerical methods are highly efficient and provide significant increase in performance over CPU-based methods when pricing PRDC swaps. An analysis of the impact of the FX volatility skew on the price of PRDC swaps is provided.
In the second part of the thesis, we develop efficient pricing algorithms for multi-asset options under the Black-Scholes-Merton framework, with strong emphasis on multi-asset American options. Our proposed pricing approach is built upon a combination of (i) a discrete penalty approach for the linear complementarity problem arising due to the free boundary and (ii) a GPU-based parallel ADI Approximate Factorization technique for the solution of the linear algebraic system arising from each penalty iteration. A timestep size selector implemented efficiently on GPUs is used to further increase the efficiency of the methods. We demonstrate the efficiency and accuracy of the proposed GPU-based parallel numerical methods by pricing American options written on three assets.
Title: | Reliable Approximate Solution of Systems of Volterra Integro-differential Equations with Time-dependent Delays |
Authors: | M. Shakourifar and W.H. Enright |
Date: | 2011 |
Pages: | 24 |
Note: | A revised version of this paper appeared in SIAM J. Sci. Comp. 33, 1134-1158. |
File: | SE_SISC.pdf |
Volterra integro-differential equations with time-dependent delay arguments (DVIDEs) can provide us with realistic models of many real-world phenomena. Delayed Lokta-Volterra predator-prey systems arise in Ecology and are well-known examples of DVIDEs first introduced by Volterra in 1928. We investigate the numerical solution of systems of DVIDEs using an adaptive stepsize selection strategy. We will present a generic variable stepsize approach for solving systems of neutral DVIDEs based on an explicit continuous Runge-Kutta method using defect error control and study the convergence of the resulting numerical method for various kind of delay arguments. We will show that the global error of the numerical solution can be effectively and reliably controlled by monitoring the size of the defect of the approximate solution and adjusting the stepsize on each step of the integration. Numerical results will be presented to demonstrate the effectiveness of this approach.
Title: | Adaptive Time-Stepping for the Strong Numerical Solution of Stochastic Differential Equations |
Authors: | Silvana Ilie, Kenneth R. Jackson and Wayne H. Enright |
Date: | January 2014. |
Pages: | 21 |
Keywords: | stochastic differential equations, adaptive time-stepping, PI control, Milstein method, strong numerical approximation, commutative noise |
Note: | Submitted to Numerical Algorithms |
File: | SDE.Adaptive.Timestepping.NA.1.pdf |
Models based on stochastic differential equations are of high interest today due to their many important practical applications. Thus the need for efficient and accurate numerical methods to approximate their solution. This paper focuses on the strong numerical solution of stochastic differential equations in Itô form, driven by multiple Wiener processes satisfying the commutativity condition. We introduce some adaptive time-stepping methods which allow arbitrary values of the stepsize. The explicit Milstein method is applied to approximate the solution of the problem and the adaptive implementations are based on a measure of the local error. Numerical tests on several models arising in applications show that our adaptive time-stepping schemes perform better than the fixed stepsize and the existing adaptive Brownian tree methods.
Title: | Calibration of Multi-Period Single-Factor Gaussian Copula Models for CDO Pricing |
Author: | Max S. Kaznady |
Date: | April, 2011 |
Pages: | 163 |
Note: | MSc. Thesis, Computer Science Dept., Univ. of Toronto, 2011. |
File: | Max_S_Kaznady_MSc_thesis_April_2011.pdf |
A Collaterized Debt Obligation (CDO) is a multi-name credit derivative, which redistributes the risk of defaults in a collection (also known as the basket or pool) of underlying assets, into fixed income securities, known as the tranches. Each tranche is associated with a certain fraction of first-to-default underlyings. Synthetic CDOs have a pool that consists of Credit Default Swaps (CDSs). If all CDSs have equal notionals, then the pool is termed homogeneous.
Single-period single-factor copula models approximate the probability of underlying defaults using a percentile to percentile transformation, and incorporate the underlying pool correlation structure for multi-name credit derivatives, such as CDOs. Currently, such models are static in time and do not calibrate consistently against market quotes. Recently Jackson, Kreinin and Zhang (JKZ) proposed a discrete-time Multi-period Single-factor Copula Model (MSCM), for which the default correlations are time-independent, allowing the model to systematically fit the market quotes. For homogeneous pools, the JKZ MSCM provides a chaining technique, which avoids expensive Monte Carlo simulation, previously used by other multi-period copula models.
However, even for homogeneous pools, the tree-based example of MSCM presented by JKZ has three drawbacks: derivatives are difficult to obtain for calibration, probabilities of the copula correlation parameter paths do not accurately represent its movements, and the model is not extremely parsimonious.
In this thesis, we develop an improved implementation of MSCM: we use an alternative multi-path parameterization of the copula correlation parameter paths and the corresponding probabilities. This allows us to calculate first-order derivatives for the MSCM in closed form for a reasonable range of parameter values, and to vary the number of parameters used by the model. We also develop and implement a practical error control heuristic for the error in the pool loss probabilities and their derivatives. We develop theoretical error bounds for the pool loss probabilities as well. We also explore a variety of optimization algorithms and demonstrate that the improved MSCM is inexpensive to calibrate. In addition, we show how MSCM calibrates to CDO data for periods before, during and after the September 2008 stock market crash.
Title: | An Efficient GPU-Based Parallel Algorithm for Pricing Multi-Asset American Options |
Authors: | Duy Minh Dang, Christina C. Christara and Kenneth R. Jackson |
Date: | March 2011 |
Pages: | 14 |
Keywords: | American Option, Multi-Asset, Penalty Method, Alternating Direction Implicit Approximate Factorization (ADI-AF), Graphics Processing Units, GPUs, Parallel Computing, Finite Difference |
Note: | A revised version of this paper (with a slightly different title) appeared in the Journal of Concurrency and Computation: Practice and Experience (JCCPE) special issue on high-performance computational finance, Vol. 24, No. 8, June 2012, pp. 849-866. |
File: | Link to the paper on the SSRN |
We develop highly-efficient parallel Partial Differential Equation (PDE) based pricing methods on Graphics Processing Units (GPUs) for multi-asset American options. Our pricing approach is built upon a combination of a discrete penalty approach for the linear complementarity problem arising due to the free boundary and a GPU-based parallel Alternating Direction Implicit Approximate Factorization technique with finite differences on uniform grids for the solution of the linear algebraic system arising from each penalty iteration. A timestep size selector implemented efficiently on GPUs is used to further increase the efficiency of the methods. We demonstrate the efficiency and accuracy of the parallel numerical methods by pricing American options written on three assets.
Title: | Pricing Multi-Asset American Options on Graphics Processing Units Using a PDE Approach |
Authors: | Duy Minh Dang, Christina C. Christara and Kenneth R. Jackson |
Date: | September 2010 |
Pages: | 8 |
Keywords: | American Option, Multi-Asset, Penalty Method, Alternating Direction Implicit Approximate Factorization (ADI-AF), Graphics Processing Units, GPUs, Parallel Computing, Finite Difference |
Note: | Proceedings of the Third Workshop on High Performance Computational Finance at SC10, The International Conference for High Performance Computing, Networking, Storage, and Analysis, sponsored by ACM/IEEE, New Orleans, USA, 13-19 November 2010. |
File: |
DOI:10.1109/WHPCF.2010.5671831 link to pdf file in ieeexplore |
We develop highly efficient parallel pricing methods on Graphics Processing Units (GPUs) for multi-asset American options via a Partial Differential Equation (PDE) approach. The linear complementarity problem arising due to the free boundary is handled by a penalty method. Finite difference methods on uniform grids are considered for the space discretization of the PDE, while classical finite differences, such as Crank-Nicolson, are used for the time discretization. The discrete nonlinear penalized equations at each timestep are solved using a penalty iteration. A GPU-based parallel Alternating Direction Implicit Approximate Factorization technique is employed for the solution of the linear algebraic system arising from each penalty iteration. We demonstrate the efficiency and accuracy of the parallel numerical methods by pricing American options written on three assets.
Title: | Energy-Based Error Control Strategies Suitable for Long MD Simulations |
Authors: | Kante Easley |
Date: | August, 2010 |
Pages: | 51 |
Note: | MSc. Thesis, Computer Science Dept., Univ. of Toronto, 2010. |
File: | Kante_Easley_Thesis.pdf |
When evaluating integration schemes used in molecular dynamics (MD) simulations, energy conservation is often cited as the primary criterion by which the integrators should be compared. As a result variable stepsize Runge-Kutta methods are often ruled out of consideration due to their characteristic energy drift. We have shown that by appropriately modifying the stepsize selection strategy in a variable-stepsize RK method it is possible for the MD practitioner to obtain substantial control over the energy drift during the course of a simulation. This ability has been previously unreported in the literature, and can be achieved without sacrificing computational efficiency under currently obtainable timescales. We present Numerical examples to confirm this.
Title: | A PDE pricing framework for cross-currency interest rate derivatives with Target Redemption features |
Authors: | Christina C. Christara, Duy Minh Dang, Kenneth R. Jackson, and Asif Lakhany |
Date: | July 2010 |
Pages: | 4 |
Keywords: | Alternating Direction Implicit, ADI, Partial Differential Equation, PDE, finite difference, Power Reverse Dual Currency swap, PRDC, Target Redemption, TARN, auxiliary state variable |
Note: | To appear in the proceedings of the International Conference of Numerical Analysis and Applied Mathematics 2010 (ICNAAM 2010), Symposium in Computational Finance, Rhodes, Greece, 19-25/09/2010. |
File: | Link to the paper on the SSRN |
We propose a general framework for efficient pricing via a partial differential equation (PDE) approach for exotic cross-currency interest rate (IR) derivatives, with strong emphasis on long-dated foreign exchange (FX) IR hybrids, namely Power Reverse Dual Currency (PRDC) swaps with a FX Target Redemption (FX-TARN) provision. The FX-TARN provision provides a cap on the FX-linked PRDC coupon amounts, and once the accumulated coupon amount reaches this cap, the underlying PRDC swap terminates. Our PDE pricing framework is based on an auxiliary state variable to keep track of the total accumulated PRDC coupon amount. Finite differences on uniform grids and the Alternating Direction Implicit (ADI) method are used for the spatial and time discretizations, respectively, of the model-dependent PDE corresponding to each discretized value of the auxiliary variable. Numerical examples illustrating the convergence properties of the numerical methods are provided.
Title: | Numerical Solution of Stochastic Models of Biochemical Kinetics |
Author: | Silvana Ilie, Wayne H. Enright and Kenneth R. Jackson |
Date: | Original draft June 2010; revised Aug. 2010 |
Pages: | 20 |
Keywords: | Chemical Master Equation, Chemical Langevin Equation, Gillespie algorithm, tau-leaping, stiff systems, multi-scale simulations, hybrid methods, stochastic biochemical modeling |
Note: | To appear in the Canadian Applied Mathematics Quarterly (CAMQ), Vol. 17, No. 3, 2009, pp. 523-554. An earlier version of this paper can be found here. |
File: | NA.Bio.Kin.Survey.rev1.pdf |
Cellular processes are typically viewed as systems of chemical reactions. Often such processes involve some species with low population numbers, for which a traditional deterministic model of classical chemical kinetics fails to accurately capture the dynamics of the system. In this case, stochastic models are needed to account for the random fluctuations observed at the level of a single cell. We survey the stochastic models of well-stirred biochemical systems and discuss important recent advances in the development of numerical methods for simulating them. Finally, we identify some key topics for future research.
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 |
Pages: | 21 |
Keywords: | Alternating Direction Implicit, ADI, Partial Differential Equation, PDE, Graphics Processing Units, GPUs, parallel computing, finite difference, multi-asset options, multi-currency swaps |
Note: | To appear in the Canadian Applied Mathematics Quarterly (CAMQ), Vol. 17, No. 4, 2009, pp. 627-660. |
File: | Link to the paper on the SSRN |
We study a parallel implementation on a Graphics Processing Unit (GPU) of Alternating Direction Implicit (ADI) time-discretization methods for solving time-dependent parabolic Partial Differential Equations (PDEs) in three spatial dimensions with mixed spatial derivatives in a variety of applications in computational finance. Finite differences on uniform grids are used for the spatial discretization of the PDEs. As examples, we apply the GPU-based parallel methods to price European rainbow and European basket options, each written on three assets. Numerical results showing the efficiency of the parallel methods are provided.
Title: | GPU Pricing of Exotic Cross-Currency Interest Rate Derivatives with a Foreign Exchange Volatility Skew Model |
Author: | Duy Minh Dang, Christina Christara and Kenneth R. Jackson |
Date: | February 2010 |
Pages: | 15 |
Keywords: | Power Reverse Dual Currency (PRDC) Swaps; Bermudan Cancelable; Partial Differential Equation (PDE); Alternating Direction Implicit (ADI), Finite Differences; Graphics Processing Units (GPUs) |
Note: | To appear in the Special Issue on "Computational Finance" of the Journal of Concurrency and Computation: Practice and Experience (JCCPE), Volume 26, Issue 9, 2014, pp.~1609--1625. |
File: | Link to the paper on the SSRN |
We present a Graphics Processing Unit (GPU) parallelization of the computation of the price of exotic cross-currency interest rate derivatives via a Partial Differential Equation (PDE) approach. In particular, we focus on the GPU-based parallel pricing of long-dated foreign exchange (FX) interest rate hybrids, namely Power Reverse Dual Currency (PRDC) swaps with Bermudan cancelable features. We consider a three-factor pricing model with FX volatility skew which results in a time-dependent parabolic PDE in three spatial dimensions. Finite difference methods on uniform grids are used 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: | Application of Generic Interpolants in the Investigation and Visualization of Approximate Solutions of PDEs on Coarse Unstructured Meshes |
Authors: | Hassan Goldani Moghaddam |
Date: | January, 2010 |
Pages: | 140 |
Note: | Ph.D. Thesis, Computer Science Dept., Univ. of Toronto, 2010. |
File: | Hassan_Goldani_Thesis.pdf |
In scientific computing, it is very common to visualize the approximate solution obtained by a numerical PDE solver by drawing surface or contour plots of all or some components of the associated approximate solutions. These plots are used to investigate the behavior of the solution and to display important properties or characteristics of the approximate solutions. In this thesis, we consider techniques for drawing such contour plots for the solution of two and three dimensional PDEs. We first present three fast contouring algorithms in two dimensions over an underlying unstructured mesh. Unlike standard contouring algorithms, our algorithms do not require a fine structured approximation. We assume that the underlying PDE solver generates approximations at some scattered data points in the domain of interest. We then generate a piecewise cubic polynomial interpolant (PCI) which approximates the solution of a PDE at off-mesh points based on the DEI (Differential Equation Interpolant) approach. The DEI approach assumes that accurate approximations to the solution and first-order derivatives exist at a set of discrete mesh points. The extra information required to uniquely define the associated piecewise polynomial is determined based on almost satisfying the PDE at a set of collocation points. In the process of generating contour plots, the PCI is used whenever we need an accurate approximation at a point inside the domain. The direct extension of the both DEI-based interpolant and the contouring algorithm to three dimensions is also investigated.
The use of the DEI-based interpolant we introduce for visualization can also be used to develop effective Adaptive Mesh Refinement (AMR) techniques and global error estimates. In particular, we introduce and investigate four AMR techniques along with a hybrid mesh refinement technique. Our interest is in investigating how well such a `generic' mesh selection strategy, based on properties of the problem alone, can perform compared with a special-purpose strategy that is designed for a specific PDE method. We also introduce an \`{a} posteriori global error estimator by introducing the solution of a companion PDE defined in terms of the associated PCI.
Title: | On Computational Methods for the Valuation of Credit Derivatives |
Author: | Wanhe Zhang |
Date: | January 2010 |
Pages: | 82 |
Note: | Ph.D. Thesis, Computer Science Dept., Univ. of Toronto, 2010. |
File: | Wanhe.Zhang.PhD.Thesis.pdf |
A credit derivative is a financial instrument whose value depends on the credit risk of an underlying asset or assets. Credit risk is the possibility that the obligor fails to honor any payment obligation. This thesis proposes four new computational methods for the valuation of credit derivatives.
Compared with synthetic collateralized debt obligations (CDOs) or basket default swaps (BDS), the value of which depends on the defaults of a prescribed underlying portfolio, a forward-starting CDO or BDS has a random underlying portfolio, as some "names" may default before the CDO or BDS starts. We develop an approach to convert a forward product to an equivalent standard one. Therefore, we avoid having to consider the default combinations in the period between the start of the forward contract and the start of the associated CDO or BDS. In addition, we propose a hybrid method combining Monte Carlo simulation with an analytical method to obtain an effective method for pricing forward-starting BDS.
Current factor copula models are static and fail to calibrate consistently against market quotes. To overcome this deficiency, we develop a novel chaining technique to build a multi-period factor copula model from several one-period factor copula models. This allows the default correlations to be time-dependent, thereby allowing the model to fit market quotes consistently. Previously developed multi-period factor copula models require multi-dimensional integration, usually computed by Monte Carlo simulation, which makes the calibration extremely time consuming. Our chaining method, on the other hand, possesses the Markov property. This allows us to compute the portfolio loss distribution of a completely homogeneous pool analytically.
The multi-period factor copula is a discrete-time dynamic model. As a first step towards developing a continuous-time dynamic model, we model the default of an underlying by the first hitting time of a Wiener process, which starts from a random initial state. We find an explicit relation between the default distribution and the initial state distribution of the Wiener process. Furthermore, conditions on the existence of such a relation are discussed. This approach allows us to match market quotes consistently.
Title: | A Survey of Modeling and Simulation of Skeletal Muscle |
Authors: | Dongwoon Lee, Michael Glueck, Azam Khan, Eugene Fiume and Ken Jackson |
Date: | December 2009 |
Pages: | 13 |
Keywords: | Computer Graphics, Three-Dimensional Graphics, Animation, Physically Based Modeling, Human Modeling, Muscle Physiology, Anatomy, Biomechanics |
Note: | Submitted to Foundations and Trends in Computer Graphics and Vision |
File: | MuscleSimulationSurvey.pdf |
Title: | Adaptive and high-order methods for valuing American options |
Authors: | Christina Christara and Duy Minh Dang |
Date: | December 2009 |
Pages: | 25 |
Keywords: | adaptive mesh selection, error equidistribution, quadratic splines, collocation, finite differences, European option, American option, penalty method |
Note: | A later version of this paper is published in the Journal of Computational Finance, Vol 14, Issue 4, 2011, pages 73-113. |
File: | ccc/americanOption.pdf |
Title: | A PDE Pricing Framework for Cross-Currency Interest Rate Derivatives |
Author: | Duy Minh Dang, Christina Christara, Kenneth R. Jackson and Asif Lakhany |
Date: | November 2009 |
Pages: | 10 |
Keywords: | Power Reverse Dual Currency Swaps, Bermudan Cancelable, Partial Differential Equation, PDE, Alternating Direction Implicit, ADI, Generalized Minimal Residual, GMRES, Fast Fourier Transform, FFT |
Note: | Accepted for publication by 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 |
Title: | Pricing of Cross-Currency Interest Rate Derivatives on Graphics Processing Units |
Author: | Duy Minh Dang |
Date: | October 2009 |
Pages: | 15 |
Keywords: | Power Reverse Dual Currency Swaps, Bermudan Cancelable, Partial Differential Equation, PDE, Alternating Direction Implicit, Finite Differences, Graphics Processing Units, Parallel Computing |
Note: | Accepted for publication in the Third International Workshop on Parallel and Distributed Computing in Finance, IEEE International Parallel & Distributed Processing Symposium, Atlanta, USA, April 2010. |
File: | Link to the paper on the SSRN |
Title: | A study of a discrete element method based granular dynamics solver |
Author: | Tina Siu Yee |
Date: | October 2009 |
Pages: | 43 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2009. |
File: | ccc/tyee-09-msc.pdf |
Keywords: | molecular dynamics, sand, powder, particle simulation, ODE solver, Runge-Kutta, friction, dampening, piling, stiction |
Title: | Quartic spline collocation for second-order boundary value problems |
Authors: | Christina C. Christara and Guohong Liu |
Date: | September 2009 |
Pages: | 8 |
Keywords: | sixth order convergence, quartic splines, second-order BVP, collocation |
Note: | Proceedings of the 9th HERCMA Conference on Computer Mathematics and Applications, Athens University of Economics and Business, September 23-26, 2009, Athens, Greece, pp. 1-8. |
File: | Please contact ccc SYMBOL_AT cs DOT toronto DOT edu |
Title: | Option Pricing Using Fourier Space Time-Stepping Framework |
Author: | Vladimir Surkov |
Date: | August 2009 |
Pages: | 134 |
Keywords: | Option pricing, Levy processes, partial integro-differential equation, Fourier transform |
Note: | Ph.D. Thesis, Computer Science Dept., Univ. of Toronto, 2009. |
File: | Link to the thesis on the SSRN |
This thesis develops a generic framework based on the Fourier transform for pricing and hedging of various options in equity, commodity, currency, and insurance markets. The pricing problem can be reduced to solving a partial integro-differential equation (PIDE). The Fourier Space Time-stepping (FST) framework developed in this thesis circumvents the problems associated with the existing finite difference methods by utilizing the Fourier transform to solve the PIDE. The FST framework-based methods are generic, highly efficient and rapidly convergent.
The Fourier transform can be applied to the pricing PIDE to obtain a linear system of ordinary differential equations that can be solved explicitly. Solving the PIDE in Fourier space allows for the integral term to be handled efficiently and avoids the asymmetrical treatment of diffusion and integral terms, common in the finite difference schemes found in the literature. For path-independent options, prices can be obtained for a range of stock prices in one iteration of the algorithm. For exotic, path-dependent options, a time-stepping methodology is developed to handle barriers, free boundaries, and exercise policies.
The thesis includes applications of the FST framework-based methods to a wide range of option pricing problems. Pricing of single- and multi-asset, European and path-dependent options under independent-increment exponential Levy stock price models, common in equity and insurance markets, can be done efficiently via the cornerstone FST method. Mean-reverting Levy spot price models, common in commodity markets, are handled by introducing a frequency transformation, which can be readily computed via scaling of the option value function. Generating stochastic volatility, to match the long-term equity options market data, and stochastic skew, observed in currency markets, is addressed by introducing a non-stationary extension of multi-dimensional Levy processes using regime-switching. Finally, codependent jumps in multi-asset models are introduced through copulas.
The FST methods are computationally efficient, running in O(MN^{d} log_{2} N) time with M time steps and N space points in each dimension on a d-dimensional grid. The methods achieve second-order convergence in space; for American options, a penalty method is used to attain second-order convergence in time. Furthermore, graphics processing units are utilized to further reduce the computational time of FST methods.
Title: | An Efficient Unified Approach for the Numerical Solution of Delay Differential Equations |
Author: | Hossein ZivariPiran and Wayne H. Enright |
Date: | September 2009 |
Pages: | 20 |
Keywords: | Delay differential equations, Discontinuous ordinary differential equations, Runge-Kutta methods, Linear multistep methods |
Note: | To appear in Numerical Algorithms DOI: 10.1007/s11075-009-9331-y Springer Link to article |
File: | hzpEnrightNA09Preprint |
In this paper we propose a new framework for designing a delay differential equation (DDE) solver which works with any supplied initial value problem (IVP) solver that is based on a standard step-by-step approach, such as Runge-Kutta or linear multi-step methods, and can provide dense output. This is done by treating a general DDE as a special example of a discontinuous IVP. Using this interpretation we develop an efficient technique to solve the resulting discontinuous IVP. We also give a more clear process for the numerical techniques used when solving the implicit equations that arise on a time step, such as when the underlying IVP solver is implicit or the delay vanishes.
The new modular design for the resulting simulator we introduce, helps to accelerate the utilization of advances in the different components of an effective numerical method. Such components include the underlying discrete formula, the interpolant for dense output, the strategy for handling discontinuities and the iteration scheme for solving any implicit equations that arise.
Title: | Efficient Simulation, Accurate Sensitivity Analysis and Reliable Parameter Estimation for Delay Differential Equations |
Author: | Hossein ZivariPiran |
Date: | September 2009 |
Pages: | 103 |
Note: | Ph.D. Thesis, Computer Science Dept., Univ. of Toronto, 2009. |
File: | hzp_phd_thesis |
Delay differential equations (DDEs) are a class of differential equations that have received considerable recent attention and been shown to model many real life problems, traditionally formulated as systems of ordinary differential equations (ODEs), more naturally and more accurately. Ideally a DDE modeling package should provide facilities for approximating the solution, performing a sensitivity analysis and estimating unknown parameters. In this thesis we propose new techniques for efficient simulation, accurate sensitivity analysis and reliable parameter estimation of DDEs.
We propose a new framework for designing a delay differential equation (DDE) solver which works with any supplied initial value problem (IVP) solver that is based on a general linear method (GLM) and can provide dense output. This is done by treating a general DDE as a special example of a discontinuous IVP. We identify a precise process for the numerical techniques used when solving the implicit equations that arise on a time step, such as when the underlying IVP solver is implicit or the delay vanishes.
We introduce an equation governing the dynamics of sensitivities for the most general system of parametric DDEs. Then, having a similar view as the simulation (DDEs as discontinuous ODEs), we introduce a formula for finding the size of jumps that appear at discontinuity points when the sensitivity equations are integrated. This leads to an algorithm which can compute sensitivities for various kind of parameters very accurately.
We also develop an algorithm for reliable parameter identification of DDEs. We propose a method for adding extra constraints to the optimization problem, changing a possibly non-smooth optimization to a smooth problem. These constraints are effectively handled using information from the simulator and the sensitivity analyzer.
Finally, we discuss the structure of our evolving modeling package DDEM. We present a process that has been used for incorporating existing codes to reduce the implementation time. We discuss the object-oriented paradigm as a way of having a manageable design with reusable and customizable components. The package is programmed in C++ and provides a user-friendly calling sequences. The numerical results are very encouraging and show the effectiveness of the techniques.
Title: | Quadratic spline collocation for one-dimensional parabolic partial differential equations |
Authors: | Christina C. Christara, Tong Chen and Dang Duy Minh |
Date: | 2009 |
Pages: | 43 |
Note: | To appear in Numerical Algorithms. |
File: | DOI
doi:10.1007/s11075-009-9317-9
Journal page: http://www.springerlink.com/content/v1g8565702373446/ Earlier version: http://www.cs.toronto.edu/pub/reports/na/ccc/parab.ps.gz |
Keywords: | quadratic splines, collocation, parabolic PDEs, Crank-Nicolson, stability, optimal order of convergence, American options. |
New methods for solving general linear parabolic partial differential equations (PDEs) in one space dimension are developed. The methods combine quadratic-spline collocation for the space discretization and classical finite differences, such as Crank-Nicolson, for the time discretization. The main computational requirements of the most efficient method are the solution of one tridiagonal linear system at each time step, while the resulting errors at the gridpoints and midpoints of the space partition are fourth order. The stability and convergence properties of some of the new methods are analyzed for a model problem. Numerical results demonstrate the stability and accuracy of the methods. Adaptive mesh techniques are introduced in the space dimension, and the resulting method is applied to the American put option pricing problem, giving very competitive results.
Title: | The Reliability/Cost Trade-off for a Class of ODE Solvers |
Authors: | W.H. Enright and L. Yan |
Date: | 2009 |
Pages: | 22 |
Note: | To appear in Numerical Algorithms DOI: 10.1007/s11075-009-9288-x Springer Link to article |
File: | enl09.pdf |
In the numerical solution of ODEs, it is now possible to develop efficient techniques that will deliver approximate solutions that are piecewise polynomials. The resulting methods can be designed so that the piecewise polynomial will satisfy a perturbed ODE with an associated defect (or residual) that is directly controlled in a consistent fashion. We will investigate the reliability/cost trade off that one faces when implementing and using such methods, when the methods are based on an underlying discrete Runge-Kutta formula.
In particular we will identify a new class of continuous Runge-Kutta methods with a very reliable defect estimator and a validity check that reflects the credibility of the estimate. We will introduce different measures of the ``reliability'' of an approximate solution that are based on the accuracy of the approximate solution; the maximum magnitude of the defect of the approximate solution; and how well the method is able to estimate the maximum magnitude of the defect of the approximate solution. We will also consider how methods can be implemented to detect and cope with special difficulties such as the effect of round-off error (on a single step) or the ability of a method to estimate the magnitude of the defect when the stepsize is large (as might happen when using a high-order method at relaxed accuracy requests).
Numerical results on a wide selection of problems will be summarized for methods of orders five, six and eight. It will be shown that a modest increase in the cost per step can lead to a significant improvement in the quality of the approximate solutions and the reliability of the method. For example, the numerical results demonstrate that, if one is willing to increase the cost per step by 50\%, then a method can deliver approximate solutions where the reported estimated maximum defect is within 1\% of its true value on 95\% of the steps.
Title: | Analytic Dynamic Factor Copula Models |
Authors: | Ken Jackson, Alex Kreinin and Wanhe Zhang |
Date: | Revised November 2011 |
Pages: | 23 |
Note: | To appear in the book Credit Portfolio Securitizations and Derivatives to be published by Wiley. Earlier versions of this paper can be found here and here. |
File: | Analytic.Dynamic.Factor.Copula.Models.2011a.pdf |
The Gaussian factor copula model is the market standard model for multi-name credit derivatives. Its main drawback is that factor copula models exhibit correlation smiles when calibrating against market tranche quotes. To overcome the calibration deficiency, we introduce a multi-period factor copula model by chaining one-period factor copula models. The correlation coefficients in our model are allowed to be time-dependent, and hence they are allowed to follow certain stochastic processes. Therefore, we can calibrate against market quotes more consistently. Usually, multi-period factor copula models require multi-dimensional integration, typically computed by Monte Carlo simulation, which makes calibration extremely time consuming. In our model, the portfolio loss of a completely homogeneous pool possesses the Markov property, thus we can compute the portfolio loss distribution analytically without multi-dimensional integration. Numerical results demonstrate the efficiency and flexibility of our model to match market quotes.
Title: | On the Blunting Method in the Verified Integration of ODEs |
Authors: | N. S. Nedialkov, K. R. Jackson and M. Neher |
Date: | July 2015; revised May 2016 |
Pages: | 20 |
Note: | Published in the Journal of Reliable Computing. The journal version can be found here. Earlier versions of this paper can be found here and here. |
File: | ned_jac_neh_blunting_2016.pdf |
Verified methods for the integration of initial value problems (IVPs) for ODEs aim at computing guaranteed error bounds for the flow of an ODE while maintaining a low level of overestimation. This paper is concerned with one of the sources of overestimation: a matrix-vector product describing a parallelepiped in phase space.
We analyze the
blunting method developed by Berz and Makino, which consists of a special
choice of the matrix in this product. For the linear model problem
u′ = Au, u(0) = u_{0} ∈ u_{0},
where u ∈ R^{m}, A ∈ R^{m x m}, m ≥ 2
and u_{0} is a given interval vector, we compare
the convergence behavior of the blunting method with that of the well-known QR
interval method.
In both methods, the amount of overestimation of the flow of the initial
set depends on the spectral radius of some well-defined matrices. We show
that under certain conditions, the spectral radii of the matrices that
describe the excess propagation in the QR method and in the blunting
method, respectively, have the same limits, so that the excess propagation
in both methods is similar.
Title: | New Interpolants for Asymptotically Correct Defect Control of BVODEs |
Authors: | W.H. Enright and P.H. Muir |
Date: | 2009 |
Pages: | 20 |
Note: | To appear in Numerical Algorithms DOI: 10.1007/s11075-009-9266-3 Springer Link to article |
File: | em09.pdf |
The defect of a continuous approximate solution to an ODE is the amount by which that approximation fails to satisfy the ODE. A number of studies have explored the use of {\it asymptotically correct} defect estimates in the numerical solution of initial value ODEs (IVODEs). By employing an appropriately constructed interpolant to an approximate discrete solution to the ODE, various researchers have shown that it is possible to obtain estimates of the local error and/or the maximum defect that are asymptotically correct on each step, as the stepsize $h \rightarrow 0$.
In this paper, we investigate the usefulness of asymptotically correct defect estimates for defect control in boundary value ODE (BVODE) codes. In the BVODE context, for a sequence of meshes which partition the problem interval, one computes a discrete numerical solution, constructs an interpolant, and estimates the maximum defect. The estimates (typically obtained by sampling the defect at a small number of points on each subinterval of the mesh) are used in a redistribution process to determine the next mesh and thus the availability of these more reliable maximum defect estimates can lead to improved meshes. As well, when such estimates are available, the code can terminate with more confidence that the defect is bounded throughout the problem domain by the user-prescribed tolerance. In this paper we employ a boot-strapping approach to derive interpolants that allow asymptotically correct defect estimates. Numerical results are included to demonstrate the validity of this approach.
Title: | Quartic spline collocation for fourth-order boundary value problems |
Authors: | Christina C. Christara, Ying Zhu and Jingrui Zhang |
Date: | September 2008 |
Pages: | 6 |
Note: | Proceedings of the 2008 Numerical Analysis conference, September 1-5, 2008, Kalamata, Greece, pp. 62-67. |
File: | Please contact ccc SYMBOL_AT cs DOT toronto DOT edu |
Title: | Multigrid and Spline Collocation Methods on Non-uniform Grids |
Author: | Guoyu Wang |
Date: | September 2008 |
Pages: | 98 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2008. |
File: | Please contact ccc SYMBOL_AT cs DOT toronto DOT edu |
Keywords: | quadratic splines, cubic splines, fourth order of convergence, collocation, multigrid, extension operator, restriction operator, preconditioning |
Title: | Randomization in the First Hitting Time Problem |
Authors: | Ken Jackson, Alex Kreinin and Wanhe Zhang |
Date: | February 2009 |
Pages: | 10 |
Note: |
A slightly revised version of this paper will appear in
Statistics and Probability Letters. An earlier version of this paper, published in Feb. 2008 and entitled "Randomization in the Default Boundary Problem", is available from http://www.defaultrisk.com/ and as First.Hitting.Time.2008.pdf. |
File: | First.Hitting.Time.2009.pdf |
In this paper we consider the following inverse problem for the first hitting time distribution: given a Wiener process with a random initial state, probability distribution, F(t), and a linear boundary, b(t) = mu t, find a distribution of the initial state such that the distribution of the first hitting time is F(t). This problem has important applications in credit risk modeling where the process represents, so-called, distance to default of an obligor, the first hitting time represents a default event and the boundary separates the healthy states of the obligor from the default state. We show that randomization of the initial state of the process makes the problem analytically tractable.
Title: | Bi-quartic Spline Collocation Methods for Fourth-order Boundary Value Problems with an Application to the Biharmonic Dirichlet Problem |
Author: | Jingrui Zhang |
Date: | January 2008 |
Pages: | 160 |
Note: | Ph.D. Thesis, Computer Science Dept., Univ. of Toronto, 2008. |
File: | Please contact ccc SYMBOL_AT cs DOT toronto DOT edu |
Keywords: | sixth order of convergence, quartic splines, collocation, fast Fourier Transform (FFT), GMRES, preconditioning |
The methods are based on quartic splines and the collocation discretization methodology with midpoints of a uniform partition as the collocation points. While the standard quartic spline method usually provides second-order approximations, the two bi-quartic spline collocation methods developed in this thesis are observed to produce approximations which are of sixth order at grid points and midpoints, and of fifth order at other points. Both are based on high order perturbations of the differential and boundary conditions operators. The one-step (extrapolated) method forces the approximations to satisfy a perturbed problem, while the three-step (deferred-correction) method proceeds in three steps and perturbs the right sides of the linear equations appropriately.
The three-step bi-quartic spline collocation method is then applied to the biharmonic Dirichlet problem and a fast solver for the resulting linear systems is developed. The fast solver consists of solving an auxiliary biharmonic problem with Dirichlet and second derivative boundary conditions along the two opposite boundaries, and a one-dimensional problem related to the two opposite boundaries. The linear system arising from the bi-quartic spline collocation discretization of the auxiliary problem is solvable by fast Fourier transforms. The problem related to the two opposite boundaries is solved by preconditioned GMRES (PGMRES). We present a detailed analysis of the performance of PGMRES by studying bounds for the eigenvalues of the preconditioned matrix. We show that the number of PGMRES iterations is independent of the grid size $n$. As a result, the complexity of the fast solver is $O(n^2\log(n))$. Numerical experiments from a variety of problems, including practical applications and problems more general than the analysis assumes, verify the accuracy of the discretization scheme and the effectiveness of the fast solver.
Title: | Fast Valuation of Forward Starting Basket Default Swaps |
Authors: | Ken Jackson, Alex Kreinin and Wanhe Zhang |
Date: | December 2007 |
Pages: | 20 |
Note: |
Also available from http://www.defaultrisk.com/ A slightly revised version of this paper is published in the International Journal of Theoretical and Applied Finance, Vol. 13, No. 2 (2010), pp. 195-209. |
File: | fwdbds.2007a.pdf |
A basket default swap (BDS) is a credit derivative with contingent payments that are triggered by a combination of default events of the reference entities. A forward-starting basket default swap (FBDS) is a BDS starting at a specified future time. Existing analytic or semi-analytic methods for pricing FBDS are time consuming due to the large number of possible default combinations before the BDS starts. This paper develops a fast approximation method for FBDS based on the conditional independence framework. The method converts the pricing of a FBDS to an equivalent BDS pricing problem and combines Monte Carlo simulation with an analytic approach to achieve an effective method. This hybrid method is a novel technique which can be viewed either as a means to accelerate the convergence of Monte Carlo simulation or as a way to estimate parameters in an analytic method that are difficult to compute directly. Numerical results demonstrate the accuracy and efficiency of the proposed hybrid method.
Title: | Analysis of the Blunting Anti-Wrapping Strategy |
Authors: | Kenneth R. Jackson, Ned S. Nedialkov and Markus Neher |
Date: | October 2007 |
Pages: | 2 |
Note: | This paper was presented at ICIAM 2007, Zurich, Switzerland, July 2007. The proceedings will be published by Proceedings in Applied Mathematics and Mechanics. |
File: | Blunting.iciam07.pdf |
Interval methods for ODEs often face two obstacles in practical computations: the dependency problem and the wrapping effect. Taylor model methods, which have been developed by Berz and his group, have recently attracted attention. By combining interval arithmetic with symbolic calculations, these methods suffer far less from the dependency problem than traditional interval methods for ODEs. By allowing nonconvex enclosure sets for the flow of a given initial value problem, Taylor model methods have also a high potential for suppressing the wrapping effect.
Makino and Berz advocate the so-called blunting method. In this paper, we analyze the blunting method (as an interval method) for a linear model ODE. We compare its convergence behavior with that of the well-known QR interval method.
Title: | Parallel Option Pricing With Fourier Space Time-Stepping Method on Graphics Processing Units |
Author: | Vladimir Surkov |
Date: | October 2007 |
Pages: | 7 |
Note: | Accepted by the First International Workshop on Parallel
and Distributed Computing in Finance 2008. The link below is to a copy of the paper on the Social Science Research Network (SSRN). |
File: | Link to the paper on the SSRN |
With the evolution of Graphics Processing Units (GPUs) into powerful and cost-efficient computing architectures, their range of application has expanded tremendously, especially in the area of computational finance. Current research in the area, however, is limited in terms of options priced and complexity of stock price models. This paper presents algorithms, based on the Fourier Space Time-stepping (FST) method, for pricing single and multi-asset European and American options with Levy underliers on a GPU. Furthermore, the single-asset pricing algorithm is parallelized to attain greater efficiency.
Title: | Numerical Methods for the Valuation of Synthetic Collateralized Debt Obligations |
Author: | Xiaofang Ma |
Date: | September 2007 |
Pages: | 82 |
Note: | Ph.D. Thesis, Computer Science Dept., Univ. of Toronto, 2007. |
File: | ma-07-phd.pdf |
A Collateralized Debt Obligation (CDO) is a credit derivative that creates fixed income securities, which are known as tranches. A CDO is called a synthetic CDO if the risky assets in the underlying pool are credit default swaps. An essential part of the valuation of a synthetic CDO tranche is how to estimate accurately and efficiently the expected value of the tranche loss function. It is known that the expected value of a function of one random variable is completely determined by the distribution of the random variable and the function itself. A standard approach to estimate the expected value of a function of one random variable is to estimate the distribution of the underlying random variable, the pool loss in our case, and then to evaluate the expected value of the given function, the tranche loss function for our problem. Following this approach, we introduce three methods for estimating the distribution of the pool loss: a stable recursive method for computing the distribution of the pool loss exactly, an improved compound Poisson approximation method and a normal power approximation method for approximating the distribution of the pool loss.
We also develop a new method that focuses on the tranche loss function directly. The tranche loss function is expressed simply in terms of two bases functions. Each of the two bases functions is a transformation of the hockey stick function h(x), where h(x) = 1 - x if 0 ≤ x < 1 and 0 if x ≥ 1. By approximating the hockey stick function by a sum of exponentials, the tranche loss function is approximated by a sum of exponentials. The main advantage of this method is that the distribution of the pool loss need not be estimated. A crucial part of this new method is the determination of the coefficients of an exponential approximation to the hockey stick function. We discuss both the numerical method for computing the exponential approximation to the hockey stick function as well as the theoretical properties of the approximation.
Performance comparisons of the four new methods developed in this thesis and other standard methods for synthetic CDO valuation are presented.
Title: | Adaptive Finite Difference Methods for Valuing American Options |
Author: | Duy Minh Dang |
Date: | September 2007 |
Pages: | 128 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2007. |
File: | Please contact ccc SYMBOL_AT cs DOT toronto DOT edu |
Keywords: | adaptive mesh selection, monitor function, finite differences, Crank-Nicolson, European option, American option, penalty method, projected SOR |
Title: | Quartic Spline Collocation Methods For Second-Order Two-Point Boundary Value ODE Problems |
Author: | Guohong Liu |
Date: | August 2007 |
Pages: | 130 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2007. |
File: | Please contact ccc SYMBOL_AT cs DOT toronto DOT edu |
Keywords: | sixth order of convergence, quartic splines, collocation, boundary conditions |
Title: | A High-Performance Method for the Biharmonic Dirichlet Problem on Rectangles |
Authors: | Christina C. Christaramand Jingrui Zhang |
Date: | July 2007 |
Pages: | 3 |
Note: | Proceedings of the 2007 International Conference On Preconditioning Techniques for Large Sparse Matrix Problems in Scientific and Industrial Applications, July 9-12, 2007 Météopole, Toulouse, France, pp. 35-37. |
File: | Please contact ccc SYMBOL_AT cs DOT toronto DOT edu |
Title: | Robust and Reliable Defect Control for Runge-Kutta Methods |
Author: | Li Yan |
Date: | July 2007 |
Pages: | 179 |
File: | LiYan-07-msc.ps.gz |
Using defect error control for a class of continuous Runge-Kutta methods for solving nonstiff problems has been studied over the last two decades. Recently a class of such methods, with an associated defect that can be reliably controlled, has been proposed \cite{Enright2007}. This approach is improved in this thesis by increasing the degree of the local interpolant by one and eliminating one of the terms contributing to the leading coefficient in the asymptotic expansion of the defect. We demonstrate this new improved robust and reliable defect approach on three continuous Runge-Kutta methods of order 5, 6, and 8. We also plot the shape of defect curves and compare them with previous work. Numerical results on the 25 test problems of DETEST at a wide range of accuracy requests are presented as well.
Title: | Option Pricing with Regime Switching Lévy Processes Using Fourier Space Time Stepping |
Authors: | Kenneth R. Jackson, Sebastian Jaimungal and Vladimir Surkov |
Date: | September 2007 (earlier draft July 2007) |
Pages: | 6 |
Note: | Appears in the proceedings of The Fourth IASTED International Conference on Financial Engineering and Applications (FEA 2007) |
File: | option_pricing_regime_switching_levy_fst_07.pdf |
Although jump-diffusion and Lévy models have been widely used in industry, the resulting pricing partial integro-differential equation poses various difficulties for valuation. Diverse finite-difference schemes for solving the problem have been introduced in the literature. Invariably, the integral and diffusive terms are treated asymmetrically, large jumps are truncated and the methods are difficult to extend to higher dimensions. We present a new efficient transform approach for regime-switching Lévy models which is applicable to a wide class of path-dependent options (such as Bermudan, barrier, and shout options) and options on multiple assets.
Title: | Efficient Contouring of PDEs on Unstructured Meshes |
Authors: | H.G. Moghaddam and W.H. Enright |
Date: | 2007 |
Pages: | 7 |
Note: | ACM Trans. on Math. Soft., 34, Article 19, 25 pages. |
File: | enm07.pdf |
We introduce three fast contouring algorithms for visualizing the solution of Partial Differential Equations based on the PCI (Pure Cubic Interpolant). The PCI is a particular piecewise bicubic polynomial interpolant defined over an unstructured mesh. Unlike standard contouring approaches, our contouring algorithms do not need a fine structured approximation and work efficiently with the original scattered data. The basic idea is to first identify the intersection points between contour curves and the sides of each triangle and then draw smooth contour segments connecting these points. We compare these contouring algorithms with the built-in Matlab procedure and with other contouring algorithms. We demonstrate that our algorithms are both more accurate and faster than the others.
Title: | Physics-Based Simulations for Fluid Mixtures |
Author: | Dongwoon Lee |
Date: | March 2007 |
Pages: | 58 |
Note: | Research Paper, Computer Science Department, University of Toronto. |
File: | DWLee.Research.Paper.2007.pdf |
We develop a physics-based particle method to simulate various fluid mixture effects. A particle method is chosen over a grid-based method because we believe a particle method is more intuitive and better suited for handling physical and chemical interactions between fluids. Since a particle carries the physical attributes of the fluid (e.g, density and velocity), transfer of these attributes between particles can be computed easily using a particle method. Also, arbitrary shapes of the fluid mixture can be modeled based on the particles' spatial configuration, which is updated throughout the simulation.
Title: | Valuation of Forward Starting Basket Default Swaps |
Author: | Wanhe Zhang |
Date: | April 2007 |
Pages: | 45 |
Note: | Research Paper, Computer Science Department, University of Toronto. |
File: | fwdbds.2007.pdf |
A basket default swap (BDS) is a credit derivative with contingent payments that are triggered by a combination of default events of the reference entities. A forward-starting BDS is a basket default swap starting at a specified future time. In this paper, we study valuation methods for a forward-starting BDS. We begin by reviewing the popular factor copula model. The widely used Monte Carlo method and associated variance reduction techniques are surveyed. The analytical solution of a recursive algorithm is developed. Conditional on a specified common factor, we explore the possible combination of defaults during the life of the forward contract; under each scenario, we evaluate the object functions; we finish the valuation by computing the expectations of these object functions. The possible combination of defaults results in a large combinatorial problem. In order to overcome the inefficiency of the method outlined above, a more applicable approximation method that omits or interpolates the unimportant scenarios is proposed. Numerical results compare the accuracy and performance of these methods and illustrate the effect of contract parameters.
Title: | Fourier Space Time-stepping for Option Pricing with Lévy Models |
Authors: | Kenneth R. Jackson, Sebastian Jaimungal and Vladimir Surkov |
Date: | March 2007 |
Pages: | 32 |
Note: | Submitted to the Journal of Computational Finance. The link below is to a copy of the paper on the Social Science Research Network (SSRN). |
File: | Link to the paper on the SSRN |
Title: | Valuation of Forward Starting CDOs |
Authors: | Ken Jackson and Wanhe Zhang |
Date: | February 2007 |
Pages: | 15 |
Note: | A slightly revised version of this paper will appear in the International Journal of Computer Mathematics |
File: | fwdcdo.2007.pdf |
Title: | Loss Distribution Evaluation for Synthetic CDOs |
Authors: | Ken Jackson, Alex Kreinin and Xiaofang Ma |
Date: | February 2007 |
Pages: | 26 |
File: | JKM.paper1.pdf |
Title: | On Convergence of Numerical Methods for Pricing Convertible Bonds |
Authors: | Dongyi (Elena) Li |
Date: | January 2007 |
Pages: | 144 |
File: | dongyi-07-msc.pdf |
In this thesis two frameworks are considered to value convertible bonds with credit risk: the TF model (Tsiveriotis and Fernandes) and the AFV model (Ayache, Forsyth and Vetzal). Both models are associated with a pair of coupled partial differential equations (PDEs) with free boundary constraints. In practice, the projected overrelaxation method and the penalty method are widely used to solve such free-boundary value problems, and the Crank-Nicolson time-stepping scheme is adopted as the underlying discretization method to achieve quadratic precision. However, due to the complexity of the PDEs in these two models and discontinuities in practice present in the boundary conditions, only linear convergence is observed in many cases. The objective of this thesis is to investigate the difficulties related to convergence and stability when solving these coupled PDEs with the Crank-Nicolson scheme, and to suggest some modifications of the basic schemes to improve stability and restore quadratic convergence.
Title: | An Exponential Approximation to the Hockey Stick Function |
Authors: | Ian Iscoe, Ken Jackson, Alex Kreinin and Xiaofang Ma |
Date: | March 2010 (replaces an earlier draft of January 2007) |
Pages: | 24 |
Note: | Submitted to the Journal of Applied Numerical Mathematics. An earlier version of this paper can be found here. |
File: | IJKM.paper2.revised.pdf |
The hockey stick (HS) function plays an important role in pricing and risk management of many financial derivatives. This paper considers approximating the HS function by a sum of exponentials. This enables the efficient computation of an approximation to the expected value of the HS function applied to a sum of conditionally independent nonnegative random variables, a task that arises in pricing many financial derivatives, CDOs in particular. The algorithm proposed by Beylkin and Monzón is used to determine the parameters of the exponential approximation to the hockey stick function. Theoretical properties of the approximation are studied and numerical results are presented.
Title: | Pricing Synthetic CDOs based on Exponential Approximations to the Payoff Function |
Authors: | Ian Iscoe, Ken Jackson, Alex Kreinin and Xiaofang Ma |
Date: | April 2011 (replaces an earlier draft of January 2007) |
Pages: | 28 |
Note: | Submitted to the Journal of Computational Finance. An earlier version of this paper can be found here. |
File: | IJKM.paper3.revised.pdf |
Correlation-dependent derivatives, such as Asset-Backed Securities (ABS) and Collateralized Debt Obligations (CDO), are common tools for offsetting credit risk. Factor models in the conditional independence framework are widely used in practice to capture the correlated default events of the underlying obligors. An essential part of these models is the accurate and efficient evaluation of the expected loss of the specified tranche, conditional on a given value of a systematic factor (or values of a set of systematic factors). Unlike other papers that focus on how to evaluate the loss distribution of the underlying pool, in this paper we focus on the tranche loss function itself. It is approximated by a sum of exponentials so that the conditional expectation can be evaluated in closed form without having to evaluate the pool loss distribution. As an example, we apply this approach to synthetic CDO pricing. Numerical results show that it is efficient.
Title: | A Scattered Data Interpolant for the Solution of Three-Dimensional PDEs |
Authors: | H.G. Moghaddam and W.H. Enright |
Date: | 2006 |
Pages: | 7 |
Note: | Proceedings of ECCOMAS CFD 2006, TU Delft, The Netherlands, September 2006. |
File: | enm06.pdf |
Using a Differential Equation Interpolant (DEI), one can accurately approximate the solution of a Partial Differential Equation (PDE) at off-mesh points. The idea is to allocate a multi-variate polynomial to each mesh element and consequently , the collection of such poly- nomials over all mesh elements will define a piecewise polynomial approximation. In this paper we will investigate such interpolants on a three-dimensional unstructured mesh. As reported elsewhere, for a tetrahedron mesh in three dimensions, tensor product tri-quadratic and pure tri-quadratic interpolants are the most appropriate candidates. We will report on the effectiveness of these alternatives on some typical PDEs.
Title: | Robust Defect Control for RK Methods Using Efficiently Computed Optimal-order Interpolants |
Authors: | W.H. Enright and W.B. Hayes |
Date: | 2006 |
Pages: | 19 |
Note: | ACM Trans. on Math. Soft., 33, pp. 1-19. |
File: | enh06.ps.gz |
The quest for reliable integration of {\it initial value problems} (IVPs) for {\it ordinary differential equations} (ODEs) is a long-standing problem in numerical analysis. At one end of the reliability spectrum are fixed stepsize methods implemented using standard floating point, where the onus lies entirely with the user to ensure the stepsize chosen is adequate for the desired accuracy. At the other end of the reliability spectrum are rigorous interval-based methods, that can provide provably correct bounds on the error of a numerical solution. This rigour comes at a price, however: interval methods are generally 2-3 orders of magnitude more expensive than fixed stepsize floating-point methods. Along the spectrum between these two extremes lie various methods of different expense that estimate and control some measure of the local errors and adjust the stepsize accordingly. In this paper, we continue previous investigations into a class of interpolants for use in Runge-Kutta methods that have a defect function whose qualitative behaviour is asymptotically independent of the problem being integrated. In particular the point, in a step, where the maximum defect occurs as $h\rightarrow 0$ is known {\it a priori}. This property allows the defect to be monitored and controlled in an efficient and robust manner even for modestly large stepsizes. Our interpolants also have a defect with the highest possible order given the constraints imposed by the order of the underlying discrete formula. We demonstrate the approach on three Runge-Kutta methods of order 5, 6, and 8, and provide Fortran and preliminary {\it Matlab} interfaces to these three new integrators. We also consider how sensitive such methods are to roundoff errors. Numerical results for four problems on a range of accuracy requests are presented.
Title: | Verifying Approximate Solutions to Differential Equations |
Authors: | W.H. Enright |
Date: | 2006 |
Pages: | 9 |
Note: | Journal of Computational and Applied Mathematics (JCAM), 185, pp. 203-311. |
File: | enr06.ps.gz |
It is now standard practice in computational science for large scale simulations to be implemented and investigated in a problem solving environment (PSE) such as MATLAB or MAPLE. In such an environment, a scientist or engineer will formulate a mathematical model, approximate its solution using an appropriate numerical method, visualize the approximate solution and verify (or validate) the quality of the approximate solution. Traditionally we have been most concerned with the development of effective numerical software for generating the approximate solution and several efficient and reliable numerical libraries are now available for use within the most widely used PSEs. On the other hand, the visualization and verification tasks have received little attention, even though each often requires as much computational effort as is involved in generating the approximate solution. In this paper we will investigate the effectiveness of a suite of tools that we have recently introduced in the MATLAB PSE to verify approximate solutions of ordinary differential equations. We will use the notion of `effectivity index', widely used by researchers in the adaptive mesh PDE community, to quantify the credibility of our verification tools. Numerical examples will be presented to illustrate the effectiveness of these tools when applied to a standard numerical method on two model test problems.
Title: | The Parallel Solution of ABD Systems Arising in Numerical Methods for BVPs for ODEs |
Author: | Richard N. Pancer |
Date: | August 2006 |
Pages: | 330 |
Note: | Ph.D. Thesis, Computer Science Dept., Univ. of Toronto, 2006. |
File: | pancer-06-phd.pdf |
Many numerical algorithms for the solution of Boundary Value Problems (BVPs) for Ordinary Differential Equations (ODEs) contain significant components that can be parallelized easily, and recent investigations have shown that substantial speedups are attainable on machines with a modest number of processors. An important step in most of these algorithms – and often the most computationally-intensive step – is the numerical solution of an Almost Block Diagonal (ABD) system of linear algebraic equations. The parallelization of this step is not so straightforward as certain characteristics of the problem make it difficult to apply standard divide-and-conquer techniques in order to arrive at a stable parallel algorithm. In the past, several algorithms have been proposed for the parallel solution of the ABD system, but many are potentially unstable or offer only limited parallelism. The proper treatment of the ABD system has proven to be a challenge in the design of parallel BVP codes.
In this thesis we present three parallel algorithms for the numerical solution of ABD systems. A parallel algorithm for this problem can be up to M / log M times faster than the fastest sequential algorithm, where the fastest sequential algorithm requires M steps. Each parallel algorithm we present attains this theoretically optimal speedup if enough processors are available, and each can be adapted for use on architectures with fewer than the required number of processors. Two of the algorithms, SLF-QR and SLF-LU, were discovered independently by us and by S.J. Wright in the 1990s. Wright presented these algorithms and analyzed their stability in two papers in the 1990s, proving SLF-QR is stable and showing that SLF-LU is stable under certain assumptions. We provide some additional insight into the stability of SLF-LU, and extend the basic algorithms to make better use of available processors during the factorization stage in order to increase parallelism in the solution stage.
The third algorithm we propose, RSCALE, is based on a notably different numerical technique: eigenvalue rescaling. RSCALE uses fewer local operations and produces less fill-in than either of the other two algorithms. In addition, RSCALE is proven to be stable for a reasonable class of problems, and has been observed to be stable for a wide class of problems through extensive numerical testing.
RSCALE is approximately 2.2 times faster than SLF-QR. The efficiency of SLF-LU is dependent on its solution strategy and its speed can vary from problem to problem, but for most problems RSCALE is approximately 1.2 times faster than SLF-LU. Moreover, we show that a variant of SLF-LU is potentially unstable on a surprising number of randomly-generated, yet well-posed, linear problems, as well as on certain nonlinear problems commonly used to test BVP codes. Since these problems pose no difficulty for either of the other two algorithms, we believe that SLF-QR, not SLF-LU, may be RSCALE's nearest competitor in terms of both speed and reliability.
Title: | Optimal Quadratic and Cubic Spline Collocation on Nonuniform Partitions |
Authors: | Christina C. Christara and Kit Sun Ng |
Date: | 2006 |
Pages: | 31 |
Note: | Appeared in Computing, 76:3-4, 2006, pp. 227-257. |
File: | DOI |
Keywords: | spline collocation; second-order two-point boundary value problem; error bounds; optimal order of convergence; spline interpolation. |
We develop optimal quadratic and cubic spline collocation methods for solving linear second-order two-point boundary value problems on non-uniform partitions. To develop optimal non-uniform partition methods, we use a mapping function from uniform to non-uniform partitions and develop expansions of the error at the non-uniform collocation points of some appropriately defined spline interpolants. The existence and uniqueness of the spline collocation approximations are shown, under some conditions. Optimal global and local orders of convergence of the spline collocation approximations and derivatives are derived, similar to those of the respective methods for uniform partitions. Numerical results on a variety of problems, including a boundary layer problem, and a non-linear problem, verify the optimal convergence of the methods, even under more relaxed conditions than those assumed by theory.
Title: | Adaptive Techniques for Spline Collocation |
Authors: | Christina C. Christara and Kit Sun Ng |
Date: | 2006 |
Pages: | 19 |
Note: | Appeared in Computing, 76:3-4, 2006, pp. 259-277. |
File: | DOI |
Keywords: | spline collocation; second-order two-point boundary value problem; error bounds; optimal order of convergence; adaptive grid; grid size estimator; error estimator; spline interpolation. |
We integrate optimal quadratic and cubic spline collocation methods for second-order two-point boundary value problems with adaptive grid techniques, and grid size and error estimators. Some adaptive grid techniques are based on the construction of a mapping function that maps uniform to non-uniform points, placed appropriately to minimize a certain norm of the error. One adaptive grid technique for cubic spline collocation is mapping-free and resembles the technique used in COLSYS (COLNEW). Numerical results on a variety of problems, including problems with boundary or interior layers, and singular perturbation problems indicate that, for most problems, the cubic spline collocation method requires less computational effort for the same error tolerance, and has equally reliable error estimators, when compared to Hermite piecewise cubic collocation. Comparison results with quadratic spline collocation are also presented.
Title: | Pricing Convertible Bonds with Dividend Protection subject to Credit Risk Using a Numerical PDE Approach |
Author: | Qingkai Mo |
Date: | April 2006 |
Pages: | 92 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2006. |
File: | ccc/moqk-06-msc.pdf |
Keywords: | convertible bond, dividend protection, credit risk, Conversion Ratio Adjustment, Dividend Pass-Thru, finite differences, Crank-Nicolson, Projected Successive Over-Relaxation (PSOR), penalty method. |
Title: | The PCI: A Scattered Data Interpolant for the Solution of Partial Differential Equations |
Authors: | H.G. Moghaddam and W.H. Enright |
Date: | 2005 |
Pages: | 10 |
Note: | Proceedings of Adaptive Modeling and Simulation 2005, (P. Diez and N.E. Wiberg eds.), pp. 243-252, Barcelona. |
File: | enhg05.ps.gz |
Most Partial Differential Equations (PDEs) that arise in practical applications do not have a closed form solution. In these cases, numerical methods can be used to approximate the solution at a discrete set of mesh points in the domain associated with the the problem definition. These approximations usually provide more accuracy at the mesh points than is available (using piecewise linear approximation, for example) at off-mesh points. Although many numerical methods have the capability of fixing the location of some or all of the mesh points, this can significantly affect the efficiency and accuracy of the computed solution. Since the best distribution of mesh points is usually non-uniform, when one wants to visualize the solution (one of the most common post-processing tasks), some extra data must be generated. Enright introduced the Differential Equation Interpolant (DEI) which approximates the solution of a PDE such that the approximations at off-mesh points have the same order of accuracy as those at mesh points.
For a two dimensional problem, the DEI is a piecewise bivariate polynomial interpolant. The number of unknown coefficients is determined by the degree and type of the interpolant (tensor product or pure) and might be greater than the number of independent linear constraints defined by the information provided by the PDE solver. Consequently, we might require extra independent constraints to uniquely determine the associated piecewise polynomial. With the DEI approach these extra constraints are based on almost satisfying the PDE at a prescribed set of 'collocation' points. The DEI generates more accurate but less smooth global approximations.
In the case that the underlying numerical method produces approximations on an unstructured mesh, the DEI can be considered to be a scattered data interpolant. Ramos and Enright introduced the Alternate Differential Equation Interpolant (ADEI), a piecewise polynomial interpolant associated with a scattered data set. Instead of choosing random points or fixed points, they developed an iterative algorithm to find suitable collocation points based on the magnitude of the coefficients of the resulting interpolant. Although the ADEI provides a relatively smooth surface, it can still suffer from the appearance of discontinuities. Moreover, finding a suitable set of collocation points may take a long time. We introduce the PCI to overcome these difficulties.
Title: | Tools for the Verification of Approximate Solutions to Differential Equations |
Authors: | W.H. Enright |
Date: | 2005 |
Pages: | 11 |
Note: | in Handbook for Scientific Computation, (B. Einarsson ed.), SIAM Press, pp. 109-119. |
File: | enr05c.ps.gz |
Scientific computation is often carried out in general purpose problem solving environments (PSEs) such as as those associated with large scale, high performance simulations. When practitioners work in these environments it is essential that they have access to state-of-the-art numerical software and tools to visualize approximate solutions produced by such software. It is equally essential (although not yet standard practice) to have available tools that can be used to verify (or validate) that the approximate solution is consistent with the true solution of the underlying mathematical model. In this presentation we will identify a suite of tools that are being developed to address this need for the case where the underlying model is a system of ODEs or PDEs. The resulting suite of tools can be used to improve the practitioners confidence in the results and allow him to quantify the inherent reliability\index{reliability}/cost trade-offs that arise. The collection of generic tools are suitable for use with any numerical method and cover a range of cost/complexity. We will assume that the approximate solution is provided as a piecewise (possibly multivariate) polynomial on an arbitrary mesh. The tools we develop perform such tasks as measuring the size of the associated defect, estimating the global error, estimating the condition number and checking that a different method generates a consistent approximate solution. The application of these tools to two typical simulations are presented.
Title: | GIA-induced secular variations in the Earth's long wavelength gravity field: Influence of 3-D viscosity variations (short communication), |
Authors: | Konstantin Latychev, Jerry X. Mitrovica, Mark E. Tamisiea, Jeroen Tromp, Christina C. Christara and Robert Moucha |
Date: | November 2005 |
Pages: | 6 |
Note: | Published in Earth and Planetary Science Letters, Vol. 240, Issue 2, Pages 322-327 |
File: | DOI link , journal page |
Keywords: | geopotential harmonics, glacial isostatic adjustment, 3-D structure |
Predictions of present day secular variations in the Earth's long wavelength geopotential driven by glacial isostatic adjustment (GIA) have previously been analyzed to infer the radial profile of mantle viscosity and to constrain ongoing cryospheric mass balance. These predictions have been based on spherically symmetric Earth models. We explore the impact of lateral variations in mantle viscosity using a new finite-volume formulation for computing the response of 3-D Maxwell viscoelastic Earth models. The geometry of the viscosity field is constrained from seismic-to-mographic images of mantle structure, while the amplitude of the lateral viscosity variations is tuned by a free parameter in the modeling. We focus on the zonal J_l harmonics for degrees l = 2, ..., 8 and demonstrate that large-scale lateral viscosity variations of two to three orders of magnitude have a modest, 5-10%, impact on predictions of J_2. In contrast, predictions of higher degree harmonics show a much greater sensitivity to lateral variation in viscosity structure. We conclude that future analyses of secular trends (for degree l > 2) estimated from ongoing (GRACE, CHAMP) satellite missions must incorporate GIA predictions based on 3-D viscoelastic Earth models.
Title: | An Efficient Algorithm Based on Quadratic Spline Collocation and Finite Difference Methods for Parabolic Partial Differential Equations |
Author: | Tong Chen |
Date: | September 2005 |
Pages: | 92 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2005. |
File: | ccc/tchen-05-msc.pdf |
Keywords: | semi-explicit semi-implicit method, optimal order of convergence, Crank-Nicolson, relaxed conditions for stability |
Title: | Pricing Convertible Bonds using Partial Differential Equations |
Author: | Lucy Xingwen Li |
Date: | September 2005 |
Pages: | 81 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2005. |
File: | lucy-05-msc.pdf |
Keywords: | Coupon payments, call, put, convert, discontinuities, implicit methods, successive over-relaxation, penalty method |
A Convertible Bond (CB) is a corporate debt security that gives the holder
the right to exchange future coupon payments and principal repayment for a
prescribed number of shares of equity. Thus, it has both an equity part
and a fixed-income part, and may contain some additional features,
such as callability and puttability.
In this paper, we study the model for valuing Convertible Bonds
with credit risk originally developed by Kostas Tsiveriotis and
Chris Fernandes (TF). The Convertible Bond is a derivative of the
stock price, and the pricing model developed by TF is based on a
free boundary value problem associated with a pair of parabolic
Partial Differential Equations (PDEs) with discontinuities at
the time points when there is a coupon payment, or when the bond
is converted, or when it is called back (purchased) by the issuer,
or when it is put (sold) to the issuer.
We explore the possible derivation of the TF model and study the
convergence of several numerical methods for solving the free
boundary value problem associated with the TF model.
In particular, we consider the Successive Over-Relaxation (SOR) iteration
and a penalty method for solving the linear complementarity problem
used to handle the free boundary.
Special emphasis is given to the effectiveness of the numerical scheme
as well to the treatment of discontinuities.
Title: | Neumaier's Method for the Solution of Initial Value Problems for Stiff Ordinary Differential Equations |
Author: | Annie Hsiao Chen Yuk |
Date: | August 2005 |
Pages: | 57 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2005. |
File: | Yuk.05.msc.pdf |
Keywords: | Stiff IVPs for ODEs, Validated Integration |
Compared with standard numerical methods for initial value problems (IVPs) for ordinary differential equations (ODEs), validated methods not only compute a numerical solution to a problem, but also generate a guaranteed bound on the global error associated with the numerical solution. There have been significant developments in the field of validated numerical methods of IVPs over the past few decades. However, none of the validated methods developed to date are suitable for stiff IVPs for ODEs.
This thesis investigates the potential of Neumaier's validated method for stiff IVPs for ODEs. We show that Neuamier's result is a special case of Dahlquist's result. Our findings show that this method has promise as an effective validated method for stiff IVPs for ODEs, for problems where there is no wrapping effect.
Title: | On Taylor Model Based Integration of ODEs |
Authors: | M. Neher, K. R. Jackson and N. S. Nedialkov |
Date: | Original Draft August 2005; revised August 2006; accepted August 2006. |
Pages: | 21 |
Note: | A slightly modified version of this paper is published in the SIAM Journal on Numerical Analysis, Vol. 45, No. 1, 2007, pp. 236-262. |
File: | Taylor_Models_ODES_2006.pdf |
Keywords: | Taylor model methods, verified integration, ODEs, IVPs |
Interval methods for verified integration of initial value problems (IVPs) for ordinary differential equations (ODEs) have been used for more than 40 years. For many classes of IVPs, these methods are able to compute guaranteed error bounds for the flow, where traditional methods provide only approximations to a solution. Overestimation, however, is a potential drawback of verified methods. For some problems, the computed error bounds become overly pessimistic, or the integration even breaks down. The dependency problem and the wrapping effect are particular sources of overestimations in interval computations.
Berz and his co-workers have developed Taylor model methods, which extend interval arithmetic with symbolic computations. The latter is an effective tool for reducing both the dependency problem and the wrapping effect. By construction, Taylor model methods appear particularly suitable for integrating nonlinear ODEs. We analyze Taylor model based integration of ODEs and compare Taylor model methods with traditional enclosure methods for IVPs for ODEs.
Title: | Glacial isostatic adjustment on 3-D Earth models: a finite-volume formulation |
Authors: | Konstantin Latychev, Jerry X. Mitrovica, Jeroen Tromp, Mark E. Tamisiea, Dimitri Komatitsch, Christina C. Christara |
Date: | May 2005 |
Pages: | 24 |
Note: | Published in Geophysical Journal International, Vol. 161, Issue 2, Pages 421-444 |
File: | DOI link , journal page |
Keywords: | crystal motions, finite volumes, glacial isostatic adjustment, 3-D viscoelastic Earth models, parallel computation, domain decomposition |
We describe and present results from a finite-volume (FV) parallel computer code for forward modelling the Maxwell viscoelastic response of a 3-D, self-gravitating, elastically compressible Earth to an arbitrary surface load. We implement a conservative, control volume discretization of the governing equations using a tetrahedral grid in Cartesian geometry and a low-order, linear interpolation. The basic starting grid honours all major radial discontinuities in the Preliminary Reference Earth Model (PREM), and the models are permitted arbitrary spatial variations in viscosity and elastic parameters. These variations may be either continuous or discontinuous at a set of grid nodes forming a 3-D surface within the (regional or global) modelling domain. In the second part of the paper, we adopt the FV methodology and a spherically symmetric Earth model to generate a suite of predictions sampling a broad class of glacial isostatic adjustment (GIA) data types (3-D crustal motions, long-wavelength gravity anomalies). These calculations, based on either a simple disc load history or a global Late Pleistocene ice load reconstruction (ICE-3G), are benchmarked against predictions generated using the traditional normal-mode approach to GIA. The detailed comparison provides a guide for future analyses (e.g. what grid resolution is required to obtain a specific accuracy?) and it indicates that discrepancies in predictions of 3-D crustal velocities less than 0.1 mm yr^{-1} are generally obtainable for global grids with ~3 x 10^6 nodes; however, grids of higher resolution are required to predict large-amplitude (>1 cm yr^{-1}) radial velocities in zones of peak post-glacial uplift (e.g. James bay) to the same level of absolute accuracy. We conclude the paper with a first application of the new formulation to a 3-D problem. Specifically, we consider the impact of mantle viscosity heterogeneity on predictions of present-day 3-D crustal motions in North America. In these tests, the lateral viscosity variation is constructed, with suitable scaling, from tomographic images of seismic S-wave heterogeneity, and it is characterized by approximately 2 orders of magnitude (peak-to-peak) lateral variations within the lower mantle and 1 order of magnitude variations in the bulk of the upper mantle (below the asthenosphere). We find that the introduction of 3-D viscosity structure has a profound impact on horizontal velocities; indeed, the magnitude of the perturbation (of order 1 mm yr^{-1}) is as large as the prediction generated from the underlying (1-D) radial reference model and it far exceeds observational uncertainties currently being obtained from space-geodetic surveying. The relative impact of lateral viscosity variations on predicted radial motions is significantly smaller.
Title: | Spline Collocation on Adaptive Grids and Non-Rectangular Domains |
Author: | Kit Sun Ng |
Date: | June 2005 |
Pages: | 187 |
Note: | Ph.D. Thesis, Computer Science Dept., Univ. of Toronto, 2005. |
File: | ngkit-05-phd.ps.gz |
Keywords: | quadratic splines, cubic splines, optimal convergence, non-uniform grid, adaptive methods, gridsize estimator, error estimator, L-shaped domains, T-shaped domains, skipped grid, moving mesh PDE. |
Optimal spline collocation is a computationally cost effective finite element method that has been relatively recently developed. Although there are many important results for optimal spline collocation, we believe the method can further be extended to solve a series of more relevant and difficult problems.
We first consider using optimal spline collocation on non-uniform partitions to solve one-dimensional linear second-order Boundary Value Problems (BVPs) that have solutions with regions of high variations, such as boundary layers. Optimal quadratic spline collocation (QSC) is extended to non-uniform partitions by using a mapping function from uniform to non-uniform partitions. Optimal cubic spline collocation (CSC) on non-uniform partitions is also examined. The formulation and analysis of the CSC methods are based, as in the QSC case, on a mapping function from uniform to non-uniform partitions. However, this method is converted into a mapping-free method and is integrated with adaptive techniques that include error and grid size estimators.
Spline collocation is then considered for solving two-dimensional linear second-order BVPs. We extend optimal spline collocation to solving BVPs on L- and T-shaped domains, which occur often in practice. The CSC method can handle general boundary conditions, which had previously never been done, but our approach requires a third step to achieve optimal convergence. We also consider two adaptive approaches for solving problems on rectangular domains whose solutions have rapid variations. As in the one-dimensional case, such problems present a significant challenge to obtain accurate approximations at minimal computational cost. We first implement a static method, which involves adding collocation points to regions where the solution has sharp variations, with CSC that uses skipped grids, which are grids that have several levels of refinement, in certain regions of the domain. An optimal skipped grid method is developed by converting the second-order BVP into a system of three first-order BVPs. We then implement a dynamic method, which involves moving the collocation points to regions where the solution has sharp variations, with CSC that incorporates the {\it Moving Mesh Partial Differential Equation} (MMPDE). The MMPDE method is further integrated into an adaptive method that can solve BVPs within a specific error tolerance.
Title: | Multigrid and Cubic Spline Collocation Methods for Advection Equations |
Author: | Zheng Zeng |
Date: | April 2005 |
Pages: | 111 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2005. |
File: | ccc/zengz-05-msc.ps.gz |
Keywords: | Shallow water equations, semi-Lagrangian methods, algebraic elimination, analytic elimination, systems of PDEs |
This thesis describes numerical methods for the solution of
advection equations in one-dimensional space. The methods are
based on combining the multigrid and cubic spline collocation (CSC)
methodologies.
Multigrid methods for cubic splines are presented.
Appropriate restriction and extension operators are developed for
cubic splines that satisfy various boundary conditions (BCs),
i.e., Dirichlet, Neumann, periodic and general BCs.
The multigrid methods are applied to CSC equations
arising from the discretization of one-dimensional
second-order differential equations.
The convergence analysis of a two-grid method integrated
with a damped Richardson relaxation scheme as smoother shows the
rate of convergence is equal to or faster than $1/2$
for a model second order equation with periodic BCs.
Multigrid methods for cubic splines are then extended to the solution of
one-dimensional shallow water equations (SWEs).
The SWEs are discretized in time with a semi-Lagrangian
semi-implicit scheme and in space with CSC methods.
At each time step, there are three different space discretization approaches:
Analytic Elimination - Discretization (AED),
Discretization -Algebraic Elimination (DAE),
and Discretization - No Elimination (DNE).
We discuss advantages and disadvantages of each, and
develop new numerical methods based on multigrid and CSC
for solving the SWEs discretized with
either the AED or the DNE approaches.
We compare our methods with Jacobi's iterative method applied
to the same problems.
Through comparison and convergence discussion, we show that our
multigrid methods for CSC are convergent and efficient.
The numerical results confirm our analysis.
Title: | DDVERK90: A User-Friendly Implementation of an Effective DDE Solver |
Author: | Hossein Zivaripiran |
Date: | April 2005 |
Pages: | 70 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2005. |
File: | hzp05.ps.gz |
Keywords: | delay differential equations, Fortran 90 solver, user-friendly solver, embedded event location, template driver. |
Delay differential equations(DDEs) provide powerful and rich mathematical models
that have been proven to be useful in simulating some real life problems.
The current publicly available DDE solvers are either very formidable to use
(especially for new users)
or designed to support only selected classes of problems.
In this thesis we have developed a new Fortran 90 DDE solver DDVERK90
that aims to combine the effectiveness of the Fortran 77 DDE solver DDVERK
and the user-friendliness of the MATLAB ODE solvers.
We have added the capability of embedded event location to
facilitate many tasks that arise in applications.
We have also introduced a new approach for classifying DDE problems
that helps users to write driver programs to solve their problem very quickly.
They do this by comparing their problem to other previously solved example problems in the same class
and using the corresponding `template' driver as a guide.
Title: | Software for Ordinary and Delay Differential Equations: Accurate Approximate Solutions are not Enough |
Authors: | W.H. Enright |
Date: | 2005 |
Pages: | 15 |
Note: | Applied Numerical Mathematics (APNUM), 56, pp. 459-471, 2006. |
File: | enr05b.pdf |
Numerical methods for both ordinary differential equations (ODEs) and delay differential equations (DDEs) are traditionally developed and assessed on the basis of how well the accuracy of the approximate solution is related to the specified error tolerance on an adaptively-chosen, discrete mesh. This may not be appropriate in numerical investigations that require visualization of an approximate solution on a continuous interval of interest (rather than at a small set of discrete points) or in investigations that require the determination of the `average' values or the `extreme' values of some solution components. In this paper we will identify modest changes in the standard error-control and stepsize-selection strategies that make it easier to develop, assess and use methods which effectively deliver approximations to differential equations (both ODEs and DDEs) that are more appropriate for these type of investigations. The required changes will typically increase the cost per step by up to 40\%, but the improvements and advantages gained will be significant. Numerical results will be presented for these modified methods applied to two example investigations (one ODE and one DDE).
Title: | Method of Lines Solutions of Boussinesq Equations |
Authors: | S. Hamdi, W.H. Enright, W.E. Schiesser and Y. Ouellet |
Date: | 2005 |
Pages: | 17 |
Note: | Journal of Computational and Applied Mathematics (JCAM), 183, pp. 327-342, 2005. |
File: | heso05.pdf |
A numerical solution procedure based on the method of lines for solving the Nwogu's one-dimensional extended Boussinesq equations is presented. The numerical scheme is accurate up to fifth order in time and fourth-order accurate in space, thus reducing all truncation errors to a level smaller than the dispersive terms retained by most extended Boussinesq models. Exact solitary wave solutions and invariants of motions recently derived by the authors are used to specify initial data for the incident solitary waves in the numerical model of Nwogu and for the verification of the associated computed solutions. The invariants of motions and several error measures are monitored in order to assess the conservative properties and the accuracy of the numerical scheme. The proposed method of lines solution procedure is general and can be easily modified to solve a wide range of Boussinesq-like equations in coastal engineering.
Title: | Exact Solutions and Conservation Laws for Generalized Quintic Regularized Long Wave Equations |
Authors: | S. Hamdi, W.H. Enright, W.E. Schiesser and J.J. Gottlieb |
Date: | 2005 |
Pages: | 10 |
Note: | Journal of Nonlinear Analysis, 63, pp. 1425-1434, 2005. Also appeared in the Proceedings of the Fourth World Congress of Nonlinear Analysis, Orlando, June 2004. |
File: | hesg05.pdf |
A new model equation is proposed for simulating unidirectional wave propagation in nonlinear media with dispersion. This model equation is a one-dimensional evolutionary partial differential equation that is obtained by coupling the generalized Korteweg-de Vries (gKdV) equation, which includes a nonlinear term of any order and cubic dispersion, with the quintic Regularized Long Wave (qRLW) equation, which includes fifth order dispersion. Exact solitary wave solutions are derived for this model equation. These analytical solutions are obtained for any order of the nonlinear term and for any given values of the coefficients of the cubic and quintic dispersive terms. Analytical expressions for three conservation laws and for three invariants of motion for solitary wave solutions of this new equation are also derived.
Title: | On the use of `arc length' and `defect' for mesh selection for differential equations |
Authors: | W.H. Enright |
Date: | 2005 |
Pages: | 7 |
Note: | to appear (2005) in Computing Letters (COLE). |
File: | enr05.ps.gz |
In this note we will consider mesh selection techniques for ordinary differential equations (ODEs) and for partial differential equations (PDEs) based on arc length and defect. We assume we are trying to approximate $y(x)$ for $x \in [0, T]$ in the ODE case and $u(x, y)$ for $(x, y) \in \Omega$ in the PDE case. The two specific areas of application that we have in mind are the numerical solution of boundary value problems (BVPs) in ODEs and the numerical solution of PDEs when a method of lines (MOL) approach is employed. The approach we develop is applicable to a wider class of PDEs including mixed-order problems and 3D problems.
Title: | An Agent-Based Simulation of Double-Auction Markets |
Author: | Ted Xiao Guo |
Date: | March 2005 |
Pages: | 75 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2005. |
File: | Guo-05.msc.pdf |
Keywords: | agent based simulation, financial markets, double-auction markets. |
Agent-based simulation has emerged as a promising research alternative to the traditional analytical approach in studying financial markets. The idea behind agent-based simulation is to constructively model a market in a bottom-up fashion using articial agents. This thesis aims to accomplish four objectives. The first objective is to present a high-quality software platform that is capable of carrying out agent-based financial market simulations. The second objective is to simulate a double-auction market with this software platform, using homogeneous zero-intelligence agents. The simulation is based on the articial market model proposed by Smith, Farmer, Gillermot and Krishnamurthy (2003). The third objective is to study the microstructure properties of the simulated double-auction market. In particular, we investigate certain statistical properties for the mid-price, bid-ask spread and limit-order book. The fourth and the last objective is to explore and evaluate trading strategies associated with the time-constrained asset liquidation problem through computational agents.
Title: | Efficient Contouring on Unstructured Meshes |
Authors: | Hassan Goldani Moghaddam |
Date: | April 2004 |
Pages: | 66 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2004. |
File: | goldani04.ps.gz |
Keywords: | interpolation, scattered data, visualization, PDE, contour, unstructured mesh |
Title: | Exact Solutions of Extended Boussinesq Equations |
Authors: | S. Hamdi, W.H. Enright, W.E. Schiesser and Y. Ouellet |
Date: | 2004 |
Pages: | 10 |
Note: | Numerical Algorithms 37, pp. 165-175. |
File: | heso04.ps.gz |
The problem addressed in this paper is the verification of numerical solutions of nonlinear dispersive wave equations such as Boussinesq-like systems of equations. A practical verification tool is to compare the numerical solution to an exact solution if available. In this work, we derive some exact solitary wave solutions and several invariants of motion for a wide range of Boussinesq-like equations using Maple software. The exact solitary wave solutions can be used to specify initial data for the incident waves in the Boussinesq numerical model and for the verification of the associated computed solution. When exact solutions are not available, four invariants of motions can be used as verification tools for the conservation properties of the numerical model.
Title: | Exact Solutions and Invariants of Motion for General types of Regularized Long Wave Equations |
Authors: | S. Hamdi, W.H. Enright, W.E. Schiesser and J.J. Gottlieb |
Date: | 2004 |
Pages: | 10 |
Note: | Mathematics and Computers in Simulation, IMACS MATCOM. |
File: | hesg04.ps.gz |
New exact solitary wave solutions are derived for general types of the regularized long wave (RLW) equation and its simpler alternative, the generalized equal width wave (EW) equation, which are evolutionary partial differential equations for the simulation of one-dimensional wave propagation in nonlinear media with dispersion processes. New exact solitary wave solutions are also derived for the generalized EW-Burgers equation, which models the propagation of nonlinear and dispersive waves with certain dissipative effects. The analytical solutions for these model equations are obtained for any order of the nonlinear terms and for any given value of the coefficients of the nonlinear, dispersive and dissipative terms. Analytical expressions for three invariants of motion for solitary wave solutions of the generalized EW equation are also devised.
Title: | A Fast Shadowing Algorithm for High Dimensional ODE Systems |
Authors: | Wayne Hayes and Kenneth R. Jackson |
Date: | November 2004 |
Pages: | 18 |
Note: | A revised version of this paper appeared in SIAM J. Sci. Comput., Vol. 29, No. 4, 2007, pp. 1738-1758. |
File: | DNA.pdf |
Keywords: | Nonlinear dynamical systems, shadowing, reliable simulation. |
Numerical solutions to chaotic dynamical systems are suspect because the chaotic nature of the problem implies that numerical errors are magnied exponentially with time. To bolster condence in the reliability of such solutions, we turn to the study of shadowing. An exact trajectory of a chaotic dynamical system lying near an approximate trajectory is called a shadow. Finding shadows of numerical trajectories of systems of ordinary dirential equations is very compute intensive, and until recently it has been infeasible to study shadows of higher-dimensional systems. This paper introduces several optimizations to previous algorithms. The optimizations collectively acheive an average speedup of almost 2 orders of magnitude with no detectable loss in ectiveness. The algorithm is tested on systems with up to 180 dimensions. Its application to large gravitational N-body integrations has already led to a deeper understanding of the reliability of galaxy simulations.
Title: | Rigorous High-Dimensional Shadowing Using Containment: The General Case |
Authors: | Carmen Young, Wayne B. Hayes, Kenneth R. Jackson |
Date: | September 2004 |
Pages: | 14 |
Note: | Submitted to the AIMS Conference on Dynamical Systems and Differential Equations |
File: | YHJ.pdf |
Keywords: | Nonlinear dynamical systems, shadowing, reliable simulation. |
A shadow is an exact solution to an iterated map that remains close to an approximate solution for a long time. An elegant geometric method for proving the existence of shadows is called containment, and it has been proven previously in two and three dimensions, and in some special cases in higher dimensions. This paper presents the general proof using tools from dirential and algebraic topology and singular homology.
Title: | Valuation of Mortgage-Backed Securities in a Distributed Environment |
Authors: | Vladimir Surkov |
Date: | March 2004 |
Pages: | 51 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2004. |
File: | Surkov04msc.pdf |
Keywords: | Valuation of Mortgage-Backed Securities, computational finance, Monte Carlo method, Quasi-Monte Carlo method, low-discrepancy sequences, Brownian bridge discretization, distributed computing environment, Microsoft .NET. |
Valuation of Mortgage-Backed Securities, regarded as integration in high-dimensional space, can be readily performed using the Monte Carlo method. The Quasi-Monte Carlo method, by utilizing low-discrepancy sequences, has been able to achieve better convergence rates for computational finance problems despite analysis suggesting that the improved convergence comes into effect only at sample sizes growing exponentially with dimension. This may be attributed to the fact that the integrands are of low effective dimension and quasi-random sequence's good equidistribution properties in low dimensions allow for the faster convergence rates to be attained at feasible sample sizes. The Brownian bridge discretization is traditionally used to reduce the effective dimension although an alternate choice of discretization can produce superior results. This paper examines the standard Brownian bridge representation and offers a reparametrization to further reduce dimension. The performance is compared both in terms of improvement in convergence and reduced effective dimensionality as computed using ANOVA decomposition. Also, porting of the valuation algorithm to a distributed environment using Microsoft .NET is presented.
Title: | Validated Numerical Bounds on the Global Error for Initial Value Problems for Stiff Ordinary Differential Equations |
Authors: | Chao Yu |
Date: | September 2004 |
Pages: | 50 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2004. |
File: | chao_yu.msc.pdf |
Keywords: | validated numerical method, initial value problem, ordinary differential equation |
There are many standard numerical methods for initial value problems (IVPs) for ordinary differential equations (ODEs). Compared with these methods, validated methods for IVPs for ODEs produce bounds that are guaranteed to contain the true solution of a problem, if the true solution exists and is unique.
The main result of this thesis is a formula to bound the global error associated with the numerical solution of a stiff IVP for an ODE. We give the complete proof of this result. Moreover, we derive Dahlquist's formula and Neumaier's formula from this formula. We also give alternative (and possibly simpler) proofs of some known related results.
Title: | Towards a high-performance method for the biharmonic Dirichlet problem |
Authors: | Jingrui Zhang |
Date: | April 2004 |
Pages: | 77 |
Note: | Research paper (equivalent to M.Sc. Thesis), Computer Science Dept., Univ. of Toronto, 2004. |
File: | ccc/jingrui-04-msc.ps.gz |
Keywords: | spline collocation; fourth-order two-point boundary value problem; optimal order of convergence; FFT; GMRES; preconditioning; eigenvalues; eigenvectors |
The discretization of the problem uses the collocation methodology and quartic splines. Quartic Spline Collocation (Quartic SC) methods that produce sixth order approximations have been recently developed for the solution of fourth order Boundary Value Problems (BVPs). In this paper, we apply an adjustment to the quartic spline basis functions so that the Quartic SC matrix has more favorable properties. Particular attention is given to the one-dimensional harmonic problem. We study the properties of two Quartic SC matrices, more specifically, one that corresponds to Dirichlet and Neumann conditions (the standard harmonic problem) and another that corresponds to Dirichlet and second derivative conditions. The latter is solvable by FFTs and used as preconditioner for the first. We develop formulae for the eigenvalues and eigenvectors of the preconditioned matrix. We show that three iterations suffice to obtain convergence for the GMRES preconditioned method. Numerical experiments verify the effectiveness of the solver and preconditioner. We also discuss the application of the preconditioned GMRES method to the two-dimensional biharmonic Dirichlet problem.
Title: | A Survey of Shadowing Methods for Numerical Solutions of Ordinary Differential Equations |
Authors: | W. B. Hayes and K. R. Jackson |
Date: | December 2003 |
Pages: | 31 |
Note: | A slightly modified version of this paper appears in the Journal Applied Numerical Mathematics, Vol. 53, 2005, pp. 299-321. |
File: | hayes-shadowing-survey.ps.gz |
A shadow is an exact solution to a set of equations that remains close to a numerical solution for a long time. Shadowing is thus a form of backward error analysis for numerical solutions to ordinary differential equations. This survey introduces the reader to shadowing with a detailed tour of shadowing algorithms and practical results obtained over the last 15 years, using the gravitational n-body problem as a running example.
Title: | Optimal Quadratic Spline Collocation on Non-Uniform Partitions |
Authors: | Christina C. Christara, and Kit Sun Ng |
Date: | June 2003, revised 2004 |
Pages: | 28 |
Note: | A revised version to appear in Computing |
File: | ccc/QSCnonuni.ps.gz |
Keywords: | spline collocation; second-order two-point boundary value problem; error bounds; optimal order of convergence; adaptive grid; grid size estimator; error estimator; spline interpolation. |
Quadratic and Cubic Spline Collocation (QSC and CSC) methods of optimal orders of convergence have been developed for the solution of linear second-order two-point Boundary Value Problems (BVPs) discretized on uniform partitions. In this paper, we extend the optimal QSC methods to non-uniform partitions, and in a companion paper we do the same for CSC. To develop optimal non-uniform partition QSC methods, we use a mapping function from uniform to non-uniform partitions and develop expansions of the error at the non-uniform collocation points of some appropriately defined spline interpolants. The existence and uniqueness of the QSC approximations are shown, under some conditions. The $j$th derivative of the QSC approximation is $O(h^{3-j})$ globally, and $O(h^{4-j})$ locally on certain points, for $j \ge 0$.
Title: | Optimal Cubic Spline Collocation on Non-Uniform Partitions |
Authors: | Christina C. Christara, and Kit Sun Ng |
Date: | June 2003, revised 2004 |
Pages: | 17 |
Note: | A revised version to appear in Computing |
File: | ccc/CSCnonuni.ps.gz |
Keywords: | spline collocation; second-order two-point boundary value problem; error bounds; optimal order of convergence; adaptive grid; grid size estimator; error estimator; spline interpolation. |
In a recent paper, non-uniform partition Quadratic Spline Collocation (QSC) methods of optimal order of convergence were developed for second-order two-point Boundary Value Problems (BVPs). In this paper, we develop optimal non-uniform Cubic Spline Collocation (CSC) methods. The formulation and analysis of the methods are based, as in the QSC case, on a mapping function from uniform to non-uniform partitions. However, the CSC method is turned into a mapping-free method. The $j$th derivative of the CSC approximation is $O(h^{4-j})$ globally, for $j \ge 0$, and $O(h^{5-j})$ locally on certain points, for $j > 0$. The non-uniform partition CSC method is integrated with adaptive grid techniques, and grid size and error estimators. One adaptive grid technique is based on the construction of a mapping function. Another is mapping-free and resembels the technique used in COLSYS (COLNEW) \cite{AMR95}. Numerical results on a variety of problems, including problems with boundary or interior layers, and singular perturbation problems indicate that, for most problems, the CSC method require less computational effort for the same error tolerance, and has equally reliable error estimators, when compared to Hermite piecewise cubic collocation (HPCC).
Title: | Fast Contouring of Solutions to Partial Differential Equations |
Authors: | E.L. Bradbury and W.H. Enright |
Date: | May 2003 |
Pages: | 25 |
Note: | This manuscript has been submitted for publication in ACM Trans. on Math. Soft., (Revised May 2003). |
File: | be03.ps.gz |
The application of Differential Equation Interpolants (DEIs) to the visualization of the solutions to Partial Differential Equations (PDEs) is investigated. In particular, we describe how a DEI can be used to generate a fine mesh approximation from a coarse mesh approximation; this fine mesh approximation can then be used by a standard contouring function to render an accurate contour plot of the surface. However, the standard approach has a time complexity equivalent to that of rendering a surface plot, $O(fm^2)$ (where $fm$ is the ratio of the width of the course mesh to the fine mesh). To address this concern three fast contouring algorithms are proposed which compute accurate contour lines directly from the DEI, and have time complexity at most $O(fm)$.
Title: | Rigorous Shadowing of Numerical Solutions of Ordinary Differential Equations by Containment |
Authors: | W. B. Hayes and K. R. Jackson |
Date: | March 2003 |
Pages: | 25 |
Note: | Accepted March 2003 for publication by the SIAM Journal on Numerical Analysis. hayes-2001-shadowing is an earlier and longer version of this paper. |
File: | hayes-2003-shadowing.ps.gz |
An exact trajectory of a dynamical system lying close to a numerical trajectory is called a shadow. We present a general-purpose method for proving the existence of finite-time shadows of numerical ODE integrations of arbitrary dimension in which some measure of hyperbolicity is present and there is either 0 or 1 expanding modes, or 0 or 1 contracting modes. Much of the rigor is provided automatically by interval arithmetic and validated ODE integration software that is freely available. The method is a generalization of a previously published containment process that was applicable only to two-dimensional maps. We extend it to handle maps of arbitrary dimension with the above restrictions, and finally to ODEs. The method involves building n-cubes around each point of the discrete numerical trajectory through which the shadow is guaranteed to pass at appropriate times. The proof consists of two steps: first, the rigorous computational verification of a simple geometric property we call the inductive containment property; and second, a simple geometric argument showing that this property implies the existence of a shadow. The computational step is almost entirely automated and easily adaptable to any ODE problem. The method allows for the rescaling of time, which is a necessary ingredient for successfully shadowing ODEs. Finally, the method is local, in the sense that it builds the shadow inductively, requiring information only from the most recent integration step, rather than more global information typical of several other methods. The method produces shadows of comparable length and distance to all currently published results.
Title: | The cost/reliability trade-off in verifying approximate solutions to differential equations |
Authors: | W.H. Enright |
Date: | February 2003 |
Pages: | 9 |
Note: | Presented at Bari International workshop on computational codes:the technological aspects of mathematics, December 2002. (Also submitted for publication in special issue of the journal of computational methods in science and engineering -- JCMSE.) |
File: | ebari03.ps.gz |
It is now standard practice in computational science for large scale simulations to be implemented and investigated in a problem solving environment (PSE) such as MATLAB or MAPLE. In such an environment, a scientist or engineer will formulate a mathematical model, approximate its solution using an appropriate numerical method, visualize the approximate solution and verify (or validate) the quality of the approximate solution. Traditionally we have been most concerned with the development of effective numerical software for generating the approximate solution and several efficient and reliable numerical libraries are now available for use within the most widely used PSEs. On the other hand, the visualization and verification tasks have received little attention, even though each often requires as much computational effort as is involved in generating the approximate solution. In this paper we will investigate the effectiveness of a suite of tools that we have recently introduced in the MATLAB PSE to verify approximate solutions of ordinary differential equations. In particular we will identify and illustrate the inherent trade-off between reliability and efficiency that arises. We will use the notion of `effectivity index', widely used by researchers in the adaptive mesh PDE community, to quantify the quality of our verification tools and illustrate the performance of these tools on a two test problems.
Title: | Exact Solutions of the Generalized Equal Width Wave Equation |
Authors: | S. Hamdi, W.H. Enright, W.E. Schiesser and J.J. Gottlieb |
Date: | March 2003 |
Pages: | 10 |
Note: | To appear in Proceedings of Workshop on Wave Phenomena in Physics and Engineering: New Models, Algorithms, and Applications, Montreal, May 2003. |
File: | hesg03.ps.gz |
The equal width wave (EW) equation is a model partial differential equation for the simulation of one-dimensional wave propagation in nonlinear media with dispersion processes. The EW-Burgers equation models the propagation of nonlinear and dispersive waves with certain dissipative effects. In this work, we derive exact solitary wave solutions for the general form of the EW equation and the generalized EW-Burgers equation with nonlinear terms of any order. We also derive analytical expressions of three invariants of motion for solitary wave solutions of the generalized EW equation.
Title: | The Design and Implementation of Usable ODE Software |
Authors: | W.H. Enright |
Date: | December 2002 |
Pages: | 13 |
Note: | Numerical Algorithms 31, 1-4, pp. 125-137, 2002. |
File: | enz02.ps.gz |
In the last decade it has become standard for students and researchers to be introduced to state-of-the-art numerical software through a problem solving environment (PSE) rather than through the use of scientific libraries callable from a high level language such as Fortran or C. In this paper we will identify the constraints and implications that this imposes on the ODE software we investigate and develop. In particular the way a `numerical solution' is displayed and viewed by a user dictates that new measures of performance and quality must be adopted. We will use the MATLAB environment and ODE software for initial value problems, boundary value problems and delay problems to illustrate the issues that arise and the progress that has been made. One of the major implications is the expectation that accurate approximations at `off-mesh' points must be provided. Traditional numerical methods for ODEs have produced approximations to the underlying solution on an associated discrete, adaptively chosen mesh. In recent years it has become common for the ODE software to also deliver approximations at off-mesh values of the independent variable. Such a feature can be extremely valuable in applications and leads to new measures of quality and performance which are more meaningful to users and more consistently interpreted and implemented in contemporary ODE software. Numerical examples of the robust and reliable behaviour of such software will be presented and the cost/reliability trade-offs that arise will be quantified.
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 |
File: | sweqsg02.ps.gz |
Note: | Submitted to Mathematics and Computers in Simulation |
Keywords: | numerical weather prediction; finite element; semi-Lagrangian semi-implicit; Rossby stability; staggered grids |
Currently in most global meteorological applications, the spectral transform method or low-order finite difference/finite element methods are used. The spectral transform method, which yields high-order approximations, requires Legendre transforms. The Legendre transforms have a computational complexity of $\mathcal O(N^3)$, where $N$ is the number of subintervals in one dimension, and thus render the spectral transform method unscalable. In this study, we present an alternative numerical method for solving the shallow water equations (SWEs) on a sphere in spherical coordinates. In this implementation, the SWEs are discretized in time using the two-level semi-Lagrangian semi-implicit method, and in space on staggered grids using the quadratic spline Galerkin method. We show that, when applied to a simplified version of the SWEs, the method yields a neutrally stable solution for the meteorologically significant Rossby waves. Moreover, we demonstrate that the Helmholtz equation arising from the discretization and solution of the SWEs should be derived algebraically rather than analytically, in order for the method to be stable with respect to the Rossby waves. These results are verified numerically using Boyd's equatorial wave equations \cite{boyd1} with initial conditions chosen to generate a soliton.
Title: | Optimal Quadratic Spline Collocation Methods for the Shallow Water Equations |
Authors: | Anita T. Layton, Christina C. Christara, and Kenneth R. Jackson |
Date: | November 2002 |
Pages: | 38 |
File: | sweqsc02.ps.gz |
Note: | Submitted to Mathematics and Computers in Simulation |
Keywords: | numerical weather prediction; finite element; semi-Lagrangian semi-implicit; Rossby stability; staggered grids |
In this study, we present numerical methods, based on the optimal quadratic spline collocation (OQSC) methods, for solving the shallow water equations (SWEs) in spherical coordinates. A quadratic spline collocation method approximates the solution of a differential problem by a quadratic spline. In the standard formulation, the quadratic spline is computed by making the residual of the differential equations zero at a set of collocation points; the resulting error is second order, while the error associated with quadratic spline interpolation is fourth order locally at certain points and third order globally. The OQSC methods generate approximations of the same order as quadratic spline interpolation. In the one-step OQSC method, the discrete differential operators are perturbed to eliminate low-order error terms, and a high-order approximation is computed using the perturbed operators. In the two-step OQSC method, a second-order approximation is generated first, using the standard formulation, and then a high-order approximation is computed in a second phase by perturbing the right sides of the equations appropriately.
In this implementation, the SWEs are discretized in time using the semi-Lagrangian semi-implicit scheme, which allows large timesteps while maintaining numerical stability, and in space using the OQSC methods. The resulting methods are efficient and yield stable and accurate representation of the meteorologically important Rossby waves. Moreover, by adopting the Arakawa C-type grid, the methods also faithfully capture the group velocity of inertia-gravity waves.
Title: | Shadowing High-Dimensional Hamiltonian Systems: the Gravitational N-Body Problem |
Author: | Wayne B. Hayes |
Date: | November 2002 |
Pages: | 4 |
Note: | Physics Review Letters, Vol. 90, No. 5, 2003. |
File: | sh-n-body-M1.ps.gz |
A shadow is an exact solution to a chaotic system of equations that remains close to a numerically computed solution for a long time. Using a variable-order, variable-timestep integrator, we numerically compute solutions to a gravitational $N$-body problem in which many particles move and interact in a fixed potential. We then search for shadows of these solutions with the longest possible duration. We find that in "softened" potentials, shadow durations are sufficiently long for significant evolution to occur. However, in unsoftened potentials, shadow durations are typically very short.
Title: | Shadowing-Based Reliability Decay in Softened N-Body Simulations |
Author: | Wayne B. Hayes |
Date: | September 2002 |
Pages: | 4 |
Note: | Astrophysical Journal Letters, Vol. 587, No. 2, Part 2, 20 April 2003, pp. L59-L62. |
File: | shadow-decay.ps.gz |
A shadow of a numerical solution to a chaotic system is an exact solution to the equations of motion that remains close to the numerical solution for a long time. In a collisionless $n$-body system, we know that particle motion is governed by the global potential rather than by inter-particle interactions. As a result, the trajectory of each individual particle in the system is independently shadowable. It is thus meaningful to measure the number of particles that have shadowable trajectories as a function of time. We find that the number of shadowable particles decays exponentially with time as $e^{-\mu t}$, and that for $\eps\in [\sim 0.2,1]$ (in units of the local mean inter-particle separation $\bar n$), there is an explicit relationship between the decay constant $\mu$, the timestep $h$ of the leapfrog integrator, the softening $\eps$, and the number of particles $N$ in the simulation. Thus, given $N$ and $\eps$, it is possible to pre-compute the timestep $h$ necessary to acheive a desired fraction of shadowable particles after a given length of simulation time. We demonstrate that a large fraction of particles remain shadowable over $\sim$100 crossing times even if particles travel up to about 1/3 of the softening length per timestep. However, a sharp decrease in the number of shadowable particles occurs if the timestep increases to allow particles to travel further than 1/3 the softening length in one timestep, or if the softening is decreased below $\sim 0.2\bar n$.
Title: | Computation of Value-at-Risk: the Fast Convolution Method, Dimension Reduction and Perturbation Theory |
Author: | Petter Wiberg |
Date: | September 2002 |
Pages: | 88 |
Note: | Ph.D. Thesis, Computer Science Dept., Univ. of Toronto, 2002 |
File: | wiberg.phd.ps.gz |
Value-at-risk is a measure of market risk for a portfolio. Market risk is the chance that the portfolio declines in value due to changes in market variables. This thesis is about the computation of value-at-risk for portfolios with derivatives and for models for returns that have a distribution with fat tails.
We introduce a new Fourier algorithm, the fast convolution method, for computing value-at-risk. The fast convolution method is different from other Fourier methods in that it does not require that the characteristic function of the portfolio returns be known explicitly. Our new method can therefore be used with more general return models. In the thesis we present experiments with three return models: the normal model, the asymmetric T model and a model using the non-parametric Parzen density estimator. We also discuss how the fast convolution method can be extended to compute the value-at-risk gradient, present a proof of convergence and illustrate the performance of the method with examples.
We develop and compare two methods for dimension reduction in the computation of value-at-risk. The goal of dimension reduction is to reduce computation time by finding a small model that captures the main dynamics of the original model. We compare the two methods for an example problem and conclude that the method based on mean square error is superior. Finally, we present an optimization example that illustrates that dimension reduction may reduce the time to compute value-at-risk while maintaining good accuracy.
We develop a perturbation theory for value-at-risk with respect to changes in the return model. By considering variational properties, we derive a first-order error bound and find the condition number of value-at-risk. We argue that the sensitivity observed in empirical studies is an inherent limitation of value-at-risk.
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 |
We revisit the optimal Quadratic Spline Collocation (QSC) methods
for single elliptic PDEs developed in \cite{Chri94}
and prove some new results, including a discrete maximum principle.
We extend QSC to
systems of two coupled linear second-order PDEs in two dimensions.
The system of PDEs is treated as one entity; no decoupling is applied.
We study the matrix properties of the linear system arising
from the discretization of systems of two PDEs by QSC.
We give sufficient conditions under which the QSC linear system
is uniquely solvable and
develop analytic formulae for the eigenvalues and eigenvectors
of the QSC matrix from $2 \times 2$ Helmholtz operators with
constant coefficients.
Numerical results demonstrate that the QSC methods are optimal,
that is, fourth order locally at certain points and third order globally,
even for some problems that do not satisfy the conditions of the analysis.
The QSC methods are compared to conventional second-order discretization
methods and are shown to produce smaller approximation errors
in the same computation time, while they achieve the same accuracy
in less time.
The formulation of the optimal two-step QSC method for
a $2 \times 2$ system of PDEs can be extended in a straightforward way to
a $n \times n$ system of PDEs.
Title: | Computation of the Probability Density Function and the Cumulative Distribution Function of the Generalized Gamma Variance Model |
Author: | Xiaofang Ma |
Date: | April 2002 |
Pages: | 77 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2002. |
File: | ma-02-msc.ps.gz |
Numerical methods for computing the probability density function (pdf) (satisfying a relative accuracy requirement) and the cumulative distribution function (cdf) (satisfying an absolute accuracy requirement) of the generalized Gamma variance model are investigated.
A hybrid method is developed to calculate the pdf. This hybrid method chooses between several basic methods to evaluate the pdf depending on the value of the risk factor and the values of the shape parameters. Extensive numerical experiments are performed to verify the robustness of the proposed hybrid method. Comparison with some existing methods for special cases suggests that this hybrid method is accurate. To improve the performance of the hybrid method, a strategy suggested by Alex Levin is adopted in the present program.
A method for computing the cdf is also developed. Numerical comparisons for several special cases are carried out to verify the correctness and the accuracy of the proposed method.
Title: | Exact Solutions of Extended Boussinesq Equations |
Authors: | S. Hamdi, W.H. Enright, Y. Ouellet and W.E. Schiesser |
Date: | December 2002 |
Pages: | 10 |
Note: | submitted to Numerical Algorithms, December 2002. |
File: | heos02.ps.gz |
The problem addressed in this paper is the verification of numerical solutions of nonlinear dispersive wave equations such as Boussinesq-like system of equations. A practical verification tool is to compare the numerical solution to an exact solution if available. In this work, we derive some exact solitary wave solutions and several invariants of motion for a wide range of Boussinesq-like equations using Maple software. The exact solitary wave solutions can be used to specify initial data for the incident waves in the Boussinesq numerical model and for the verification of the associated computed solution. When exact solutions are not available, four invariants of motions can be used as verification tools for the conservation properties of the numerical model.
Title: | Interpolation of Surfaces over Scattered Data |
Authors: | G.A. Ramos and W.H. Enright |
Date: | September 2001 |
Pages: | 7 |
Note: | Proceedings of IASTED, (M.H. Hamza, ed.), pp.219-224, Marbella Spain, September 2001, ACTA Press. |
File: | er01.ps.gz |
We investigate the performance of DEI, an approach [2] that computes off-mesh approximations of PDE solutions, and that can also be used as a technique for scattered data interpolation and surface representation. For the general case of unstructured meshes, we found it necessary to modify the original DEI. The resulting method, ADEI, adjusts the parameter of the interpolant, obtaining better performance. Finally we measure ADEI's performance using different sets of scattered data and compare its performance against two methods from the collection of ACM algorithms.
Title: | Hierarchichal Representation and Computation of Approximate Solutions in Scientific Simulations |
Authors: | W.H. Enright |
Date: | 2001 |
Pages: | 17 |
Note: | in `The Architecture of Scientific Software', (R. Boisvert and P. Tang eds.), pp. 301-316, Kluwer, Boston. |
File: | ebt01.ps.gz |
In many applications involving large scale scientific computing a massive amount of data is generated in a typical simulation. Such simulations generally require the numerical solution of systems of differential equations where the results are often generated remotely using special high-performance software and computer systems visualization tools. The visualization packages are usually run on local workstations and make use of colour, lighting, texture, sound and animation to highlight and reveal interesting characteristics or features of the approximate solution. This `interactive' viewing of the data is the `rendering' stage of the modeling process and it can be very selective and local in the sense that only a subset of the variables are rendered and then only in regions where something interesting is happening. The result is that, in many simulations, large amounts of data must be stored and shared in a distributed environment while only a very small fraction of this data will ever be viewed by anyone. In this paper we propose an approach, using a hierarchichal representation of the approximate solution, which will avoid the generation and storage of data that is not required.
Title: | Rigorous Shadowing of Numerical Solutions of Ordinary Differential Equations by Containment |
Authors: | W. B. Hayes and K. R. Jackson |
Date: | December 2001 |
Pages: | 41 |
Note: | hayes-2003-shadowing is a shorter, revised version of this paper. |
File: | hayes-2001-shadowing.ps.gz |
An exact trajectory of a dynamical system lying close to a numerical trajectory is called a shadow. We present a general-purpose method for proving the existence of finite-time shadows of numerical ODE integrations of arbitrary dimension in which some measure of hyperbolicity is present and there is either 0 or 1 expanding modes, or 0 or 1 contracting modes. Much of the rigor is provided automatically by interval arithmetic and validated ODE integration software that is freely available. The method is a generalization of a previously published containment process that was applicable only to two-dimensional maps. We extend it to handle maps of arbitrary dimension with the above restrictions, and finally to ODEs. The method involves building n-cubes around each point of the discrete numerical trajectory through which the shadow is guaranteed to pass at appropriate times. The proof consists of two steps: first, the rigorous computational verification of a simple geometric property we call the inductive containment property; and second, a simple geometric argument showing that this property implies the existence of a shadow. The computational step is almost entirely automated and easily adaptable to any ODE problem. The method allows for the rescaling of time, which is a necessary ingredient for successfully shadowing ODEs. Finally, the method is local, in the sense that it builds the shadow inductively, requiring information only from the most recent integration step, rather than more global information typical of several other methods. The method produces shadows of comparable length and distance to all currently published results.
Title: | Dimension Reduction in the Computation of Value-at-Risk |
Authors: | Claudio Albanese, Ken Jackson and Petter Wiberg |
Date: | November 2001 |
Pages: | 19 |
Note: | Accepted May 2002 for publication in the Journal of Risk Finance |
File: | dimred.ps.gz |
Value-at-risk models can have many dimensions. We present two new algorithms for dimension reduction in value-at-risk algorithms with Delta-Gamma approximations. In the first method, we compute a reduced portfolio with a small mean square error for the residual and, in the second method, we use low rank approximations to find a reduced portfolio. The paper concludes with an example, estimating value-at-risk and hedging an option portfolio, that demonstrates that dimension reduction leads to large savings in computational time without sacrificing accuracy.
Title: | The Design and Implementation of an Object-Oriented Validated ODE Solver |
Authors: | N. S. Nedialkov and K. R. Jackson |
Date: | Original draft October 2001; Revised June 2002. |
Pages: | 32 |
File: | ned.software.01.ps.gz |
Validated methods for initial value problems (IVPs) for ordinary differential equations (ODEs) produce bounds that are guaranteed to contain the true solution of a problem. These methods use interval techniques to enclose rounding errors, truncation errors, and uncertainties in the model.
We describe the design and implementation of VNODE, Validated Numerical ODE, which is an object-oriented solver for computing rigorous bounds on the solution of an IVP for an ODE. We have designed VNODE as a set of C++ classes, which encapsulate data structures and algorithms employed in validated methods for IVPs for ODEs.
In VNODE, new solvers can be constructed by choosing methods from sets of methods. In addition, part of a solver can be replaced by another one performing the same task. This enables the user to isolate and compare methods implementing the same part of a solver. Moreover, a new method can be added to VNODE by encapsulating it in a class and deriving this class from an existing one.
We illustrate how to integrate with VNODE, how to construct new solvers, and how to extend this package with new methods.
Title: | Quartic-Spline Collocation Methods for Fourth-Order Two-Point Boundary Value Problems |
Author: | Ying Zhu |
Date: | April 2001 |
Pages: | 73 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 2001. |
File: | ccc/yz-01-msc.ps.gz |
This thesis presents numerical methods for the solution of
general linear fourth-order boundary value problems in one dimension.
The methods are based on quartic splines, that is, piecewise quartic
polynomials with C^3 continuity, and the collocation discretization
methodology with the midpoints of a uniform partition being the
collocation points.
The standard quartic-spline collocation method is shown
to be second order. Two sixth-order quartic-spline collocation methods
are developed and analyzed. They are both based on a high order
perturbation of the differential equation and boundary conditions
operators. The one-step method forces the approximation to satisfy
a perturbed problem, while the three-step method proceeds in three steps
and perturbs the right sides of the equations appropriately.
The error analysis follows the Green's function approach and shows
that both methods exhibit optimal order of convergence, that is,
they are locally sixth order on the gridpoints and midpoints
of the uniform partition, and fifth order globally.
The properties of the matrices arising from a restricted class
of problems are studied. Analytic formulae for the eigenvalues
and eigenvectors are developed, and related to those arising
from quadratic-spline collocation matrices.
Numerical results verify the orders of convergence predicted by the
analysis.
Title: | High-Order Spatial Discretization Methods for the Shallow Water Equations |
Author: | Anita W. Tam |
Date: | February 2001 |
Pages: | 160 |
Note: | Ph.D. Thesis, Computer Science Dept., Univ. of Toronto, 2001. |
File: | tam-01-phd.ps.gz |
We present new numerical methods for the shallow water equations on a sphere in spherical coordinates. In our implementation, the equations are discretized in time with the two-level semi-Lagrangian semi-implicit (SLSI) method, and in space on a staggered grid with the quadratic spline Galerkin (QSG) and the optimal quadratic spline collocation (OQSC) methods. When discretized on a uniform spatial grid, the solutions are shown through numerical experiments to be fourth-order in space locally at the nodes and midpoints of the spatial grids, and third-order globally.
We also show that, when applied to a simplified version of the shallow water equations, each of our algorithms yields a neutrally stable solution for the meteorologically significant Rossby waves. Moreover, we demonstrate that the Helmholtz equation associated with the shallow water equations should be derived algebraically rather than analytically in order for the algorithms to be stable with respect to the Rossby waves. These results are verified numerically using Boyd's equatorial wave equations with initial conditions to generate a soliton.
We then analyze the performance of our methods on various staggered grids --- the A-, B-, and C-grids. From an eigenvalue analysis of our simplified version of the shallow water equations, we conclude that, when discretized on the C-grid, our algorithms faithfully capture the group velocity of inertia-gravity waves. Our analysis suggests that neither the A- nor B-grids will produce such good results. Our theoretical results are supported by numerical experiments, in which we discretize Boyd's equatorial wave equations using different staggered grids and set the initial conditions for the problem to generate gravitation modes instead of a soliton. With both the A- and B-grids, some waves are observed to travel in the wrong direction, whereas, with the C-grid, gravity waves of all wavelengths propagate in the correct direction.
Title: | Hedging with Value-At-Risk |
Author: | Claudio Albanese, Ken Jackson and Petter Wiberg |
Date: | January 2001 |
Pages: | 17 |
Note: | |
File: | vargrad.ps.gz |
Methods to computing value-at-risk gradients with respect to portfolio positions have many applications. They include calculation of capital/reward efficient frontiers, hedging of derivative portfolios and optimal replication. We present a new algorithm for computing value-at-risk and its gradients. If the return can be decomposed as a sum of independent portfolio marginals, the pay-off distribution can be computed with multiple convolutions. This principal-component-type decomposition is also useful for calibrating fat-tailed distributions. We conclude with two applications of hedging with value-at-risk.
Title: | Some Recent Advances in Validated Methods for IVPs for ODEs |
Author: | N. S. Nedialkov and K. R. Jackson |
Date: | January 2001 |
Pages: | 7 |
Note: | To appear in the Proceedings of the International Minisymposium on Mathematical Modelling and Scientific Computations, 8-11 April 2001, Borovets, Bulgaria. |
File: | ned.bg.01.ps.gz |
We present an overview of interval Taylor series (ITS) methods for IVPs for ODEs and discuss some recent advances in the theory of validated methods for IVPs for ODEs, such as, an interval Hermite-Obreschkoff (IHO) scheme, the stability of ITS and IHO methods, and a new perspective on the wrapping effect.
Title: | Rigorous Shadowing of Numerical Solutions of Ordinary Differential Equations by Containment |
Author: | Wayne B. Hayes |
Date: | January 2001 |
Pages: | 87 |
Note: | Ph.D. Thesis, Computer Science Dept., Univ. of Toronto, 2001. |
File: | hayes-01-phd.ps.gz |
An exact trajectory of a dynamical system lying close to a numerical trajectory is called a shadow. We present a general-purpose method for proving the existence of finite-time shadows of numerical ODE integrations of arbitrary dimension in which some measure of hyperbolicity is present and there is either 0 or 1 expanding modes, or 0 or 1 contracting modes. Much of the rigor is provided automatically by interval arithmetic and validated ODE integration software that is freely available. The method is a generalization of a previously published containment process that was applicable only to two-dimensional maps. We extend it to handle maps of arbitrary dimension with the above restrictions, and finally to ODEs. The method involves building n-cubes around each point of the discrete numerical trajectory through which the shadow is guaranteed to pass at appropriate times. The proof consists of two steps: first, the rigorous computational verification of an inductive containment property; and second, a simple geometric argument showing that this property implies the existence of a shadow. The computational step is almost entirely automated and easily adaptable to any ODE problem. The method allows for the rescaling of time, which is a necessary ingredient for successfully shadowing ODEs. Finally, the method is local, in the sense that it builds the shadow inductively, requiring information only from the most recent integration step, rather than more global information typical of several other methods. The method produces shadows of comparable length and distance to all currently published results.
Title: | Superconvergent Interpolants for Collocation Methods applied to mixed order BVODEs |
Authors: | W.H. Enright and R. Sivasothinathan |
Date: | 2000 |
Pages: | 34 |
Note: | A revised version of this paper appeared in ACM Trans. on Math. Soft.,26, 3, pp.323-35. |
File: | er00.ps.gz |
Continuous approximations to boundary value problems in ordinary differential equations (BVODEs), constructed using collocation at Gauss points, are more accurate at the mesh points than at off-mesh points. From these approximations, it is possible to construct improved continuous approximations by extending the high accuracy that is available at the mesh points to off-mesh points. One possibility is the {Bootstrap} approach, which improves the accuracy of the approximate solution at the off-mesh points in a sequence of steps until the accuracy at the mesh points and off-mesh points is consistent. A Bootstrap approach for systems of mixed order BVODEs is developed to improve approximate solutions produced by {COLNEW}, a Gauss collocation based software package. An implementation of this approach is discussed and numerical results presented which confirm that the improved approximations satisfy the predicted error bounds and are relatively inexpensive to construct.
Title: | Continuous Numerical Methods for ODEs with Defect Control |
Authors: | W.H. Enright |
Date: | 2000 |
Pages: | 16 |
Note: | Journal of Computational and Applied Mathematics, 125 (2000), pp. 159-170. |
File: | edc01.ps.gz |
Over the last decade several general purpose numerical methods for ordinary differential equations (ODEs) have been developed which generate a continuous piecewise polynomial approximation that is defined for all values of the independent variable in the range of interest. For such methods it is possible to introduce measures of the quality of the approximate solution based on how well the piecewise polynomial satisfies the ODE. This leads naturally to the selection strategies in order to control the magnitude of the associated defect can be very effective and such methods are now being widely used. In this paper we will review the advantages of this class of numerical methods and present examples of how they can be effectively applied. We will focus on numerical methods for initial value problems (IVPs) and boundary value problems (BVPs) where most of the developments have been introduced but we will also discuss the implications and related developments for other classes of ODEs such as delay differential equations (DDEs) and differential algebraic equations (DAEs).
Title: | Fast Fourier Transform Solvers and Preconditioners for Quadratic Spline Collocation |
Authors: | Christina C. Christara and Kit Sun Ng |
Date: | November 2000 (rev. October 2001, February 2002) |
Pages: | 33 |
Note: | Depart. of Computer Science, University of Toronto, Technical Report 317/2000. A revised version to appear in BIT, 42:4 (December 2002). |
File: | ccc/precond.ps.gz |
DOI: | 10.1023/A:1021944218806 |
Journal: | PDF file (fulltext), page |
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, along at least one direction of a rectangular domain. General variable coefficient PDEs are handled by preconditioned iterative solvers. The preconditioner is the QSC matrix arising from a constant coefficient PDE. The convergence analysis of the preconditioner is presented. It is shown that, under certain conditions, the convergence rate is independent of the gridsize. The preconditioner is solved by FFT techniques, and integrated with one-step or acceleration methods, giving rise to asymptotically almost optimal linear solvers, with complexity O(N^2 log N). Numerical experiments verify the effectiveness of the solvers and preconditioners, even on problems more general than the analysis assumes. The development and analysis of FFT solvers and preconditioners is extended to QSC equations corresponding to systems of elliptic PDEs.
Title: | A New Perspective on the Wrapping Effect in Interval Methods for Initial Value Problems for Ordinary Differential Equations |
Authors: | Nedialko S. Nedialkov and Kenneth R. Jackson |
Date: | November 2000 |
Pages: | 46 |
Note: | A slightly modified version of this paper is published in the book "Perspectives on Enclosure Methods", Ulrich Kulisch, Rudolf Lohner and Axel Facius, editors, Springer-Verlag, 2001, pp. 219-264. |
File: | ned.scan00.ps.gz |
The problem of reducing the wrapping effect in interval methods for initial value problems for ordinary differential equations has usually been studied from a geometric point of view. We develop a new perspective on this problem by linking the wrapping effect to the stability of the interval method. Thus, reducing the wrapping effect is related to finding a more stable scheme for advancing the solution. This allows us to exploit eigenvalue techniques and to avoid the complicated geometric arguments used previously.
We study the stability of several anti-wrapping schemes, including Lohner's QR-factorization method, which we show can be viewed as a simultaneous iteration for computing the eigenvalues of an associated matrix. Using this connection, we explain how Lohner's method improves the stability of an interval method and show that, for a large class of problems, its global error is not much bigger than that of the corresponding point method.
Title: | Some Recent Advances in Validated Methods for IVPs for ODEs |
Authors: | Kenneth R. Jackson and Nedialko S. Nedialkov |
Date: | Revised August 2001 (original report November 2000) |
Pages: | 23 |
Note: | Accepted December 2001 for publication in the Journal of Applied Numerical Mathematics. |
File: | numdiff9a.ps.gz |
Compared to standard numerical methods for initial value problems (IVPs) for ordinary differential equations (ODEs), validated methods (often called interval methods) for IVPs for ODEs have two important advantages: if they return a solution to a problem, then (1) the problem is guaranteed to have a unique solution, and (2) an enclosure of the true solution is produced.
We present a brief overview of interval Taylor series (ITS) methods for IVPs for ODEs and discuss some recent advances in the theory of validated methods for IVPs for ODEs. In particular, we discuss an interval Hermite-Obreschkoff (IHO) scheme for computing rigorous bounds on the solution of an IVP for an ODE, the stability of ITS and IHO methods, and a new perspective on the wrapping effect, where we interpret the problem of reducing the wrapping effect as one of finding a more stable scheme for advancing the solution.
Title: | PMIRKDC: a Parallel Mono-Implicit Runge-Kutta Code with Defect Control for Boundary Value ODEs |
Authors: | P.H. Muir, R.N. Pancer, K.R. Jackson |
Date: | Original Draft September 2000; Revised January 2003. |
Pages: | 34 |
Note: | A slightly modified version of this paper appeared in
Parallel Computing, Volume 29, Issue 6, June 2003,
Pages 711-741. The following report is a longer version of this paper. |
File: | parallel-mirkdc.ps.gz |
We describe parallel software, PMIRKDC, for solving boundary value ordinary differential equations (BVODEs). This software is based on the package, MIRKDC, which employs mono-implicit Runge-Kutta schemes within a defect control algorithm. The primary computational costs involve the treatment of large, almost block diagonal (ABD) linear systems. The most significant feature of PMIRKDC is the replacement of sequential ABD software, COLROW, with new parallel ABD software, RSCALE, based on a parallel block eigenvalue rescaling algorithm. Other modifications involve parallelization of the setup of the ABD systems and solution interpolants, and defect estimation. Numerical results show almost linear speedups.
Title: | Runge-Kutta Software for the Parallel Solution of Boundary Value ODEs |
Authors: | P.H. Muir, R.N. Pancer, K.R. Jackson |
Date: | August 2000 |
Pages: | 46 |
Note: | The previous paper is a shorter version of this technical report. |
File: | parallel-mirkdc-tr.ps.gz |
In this paper we describe the development of parallel software for the numerical solution of boundary value ordinary differential equations (BVODEs). The software, implemented on two shared memory, parallel architectures, is based on a modification of the MIRKDC package, which employs discrete and continuous mono-implicit Runge-Kutta schemes within a defect control algorithm. The primary computational costs are associated with the almost block diagonal (ABD) linear systems representing the Newton matrices arising from the iterative solution of the nonlinear algebraic systems which result from the discretization of the ODEs. The most significant modification featured in the parallel version of the code is the replacement of the sequential ABD linear system software, COLROW, which employs alternating row and column elimination, with new parallel ABD linear system software, RSCALE, which is based on a recently developed parallel block eigenvalue rescaling algorithm. Other modifications are associated with the parallelization of the setup of the ABD systems, the setup of the approximate solution interpolants, and the estimation of the defect. The numerical results show that nearly optimal speedups can be obtained, and that substantial speedups in overall solution time are achieved, compared with the sequential version of MIRKDC.
Title: | Quadratic Spline Collocation Methods for Systems of Elliptic PDEs |
Author: | Kit Sun Ng |
Date: | April 2000 |
Pages: | 67 |
Note: | M.Sc. Thesis, Computer Science Departmnent, University of Toronto, 2000. |
File: | ccc/ngkit-00-msc.ps.gz |
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: | An Effective High-Order Interval Method for Validating Existence and Uniqueness of the Solution of an IVP for an ODE |
Authors: | N. S. Nedialkov, K. R. Jackson and J. D. Pryce |
Date: | October 2000. (This is a revision of our original draft of December 1999.) |
Pages: | 17 |
Note: | Reliable Computing, Vol. 7, Issue 6, Nov. 2001, pp. 449-465. |
File: | stepsize.99a.ps.gz |
Validated methods for initial value problems for ordinary differential equations produce bounds that are guaranteed to contain the true solution of a problem. When computing such bounds, these methods verify that a unique solution to the problem exists in the interval of integration and compute a priori bounds for the solution in this interval. A major difficulty in this verification phase is how to take as large a stepsize as possible, subject to some tolerance requirement. We propose a high-order enclosure method for proving existence and uniqueness of the solution and computing a priori bounds.
Title: | Computing Rigorous Bounds on the Solution of an Initial Value Problem for an Ordinary Differential Equation |
Authors: | Nedialko (Ned) Stoyanov Nedialkov |
Date: | February 1999 |
Pages: | 149 |
Note: | Ph.D. Thesis, Computer Science Dept., Univ. of Toronto, 1999. |
File: | ned-99-phd.ps.gz |
Compared to standard numerical methods for initial value problems (IVPs) for ordinary differential equations (ODEs), validated (also called interval) methods for IVPs for ODEs have two important advantages: if they return a solution to a problem, then (1) the problem is guaranteed to have a unique solution, and (2) an enclosure of the true solution is produced.
To date, the only effective approach for computing guaranteed enclosures of the solution of an IVP for an ODE has been interval methods based on Taylor series. This thesis derives a new approach, an interval Hermite-Obreschkoff (IHO) method, for computing such enclosures.
Compared to interval Taylor series (ITS) methods, for the same order and stepsize, our IHO scheme has a smaller truncation error and better stability. As a result, the IHO method allows larger stepsizes than the corresponding ITS methods, thus saving computation time. In addition, since fewer Taylor coefficients are required by IHO than ITS methods, the IHO method performs better than the ITS methods when the function for computing the right side contains many terms.
The stability properties of the ITS and IHO methods are investigated. We show as an important by-product of this analysis that the stability of an interval method is determined not only by the stability function of the underlying formula, as in a standard method for an IVP for an ODE, but also by the associated formula for the truncation error.
This thesis also proposes a Taylor series method for validating existence and uniqueness of the solution, a simple stepsize control, and a program structure appropriate for a large class of validated ODE solvers.
Title: | An Efficient Transposition Algorithm for Distributed Memory Computers |
Authors: | C. C. Christara, X. Ding and K. R. Jackson |
Date: | October 1999. (This is a revision of our drafts of January and April 1999.) |
Pages: | 19 |
Note: | Accepted for publication in the proceedings of the 13th Annual International Symposium on High Performance Computing Systems and Applications (HPCS'99), Queen's University, Kingston, Ontario, Canada, 13-16 June 1999. Appeared as Chapter 38 in High Performance Computing Systems and Applications 1999, A. Pollard, D. J. K. Mewhort and D. F. Weaver eds. pp. 349-368, Kluwer Academic Publishers. |
File: | transpose.99b.ps.gz |
Title: | ODE Software that Computes Guaranteed Bounds on the Solution |
Authors: | N. S. Nedialkov and K. R. Jackson |
Date: | Original draft December 1998; Revised March 1999. |
Pages: | 31 |
Note: | Accepted as a chapter in the book
"Advances in Software Tools for Scientific Computing",
edited by H. P. Langtangen, A. M. Bruaset and E. Quak,
Springer-Verlag, 1999. The title of the original draft was "Software Issues in Validated ODE Solving". |
File: | SciTools.98.ps.gz |
Title: | An Interval Hermite-Obreschkoff Method for Computing Rigorous Bounds on the Solution of an Initial Value Problem for an Ordinary Differential Equation |
Authors: | N. S. Nedialkov and K. R. Jackson |
Date: | February 1999 (original draft October 1998) |
Pages: | 22 |
Note: | Appeared in Reliable Computing, 5 (1999), pp. 289-310 and as a chapter in the book Developments in Reliable Computing, T. Csendes, editor, pages 289-310, Kluwer, Dordrecht, Netherlands, 1999. |
File: | ned.scan98.ps.gz |
The stability properties of the ITS and IHO methods are investigated. We show as an important by-product of this analysis that the stability of an interval method is determined not only by the stability function of the underlying formula, as in a standard method for an IVP for an ODE, but also by the associated formula for the truncation error.
Title: | Numerical Solution of the Shallow-Water Equations on Distributed Memory Systems |
Author: | Xiaoliang (Lloyd) Ding |
Date: | October 1998 |
Pages: | 73 |
Note: | M.Sc. Thesis, Computer Science Departmnent, University of Toronto, 1998. |
File: | ding.98.msc.ps.gz |
The data transposition requires an all-to-all type communication. The Direct Exchange algorithm, commonly used for this task, is inefficient if the number of processors is large. We investigate a series of more sophisticated techniques: the Ring Exchange, Mesh Exchange and Cube Exchange algorithms. We compare the performance of these techniques with that of the Message Passing Interface collective communication routine. We perform an isoefficiency analysis to measure the scalability of the parallel system. The numerical experiments are carried out on a Cray T3E with 512 processors and on an ethernet-connected cluster of 36 workstations.
Both the analysis and our numerical results indicate that the Cube Exchange and Mesh Exchange algorithms perform better than the other methods. However, even with the more sophisticated techniques, the parallel algorithm for the shallow-water equations remains essentially unscalable.
Title: | A Numerical Study of One-Factor Interest Rate Models |
Author: | Zhong Ge |
Date: | January 1998 |
Pages: | 95 |
Note: | M.Sc. Thesis, Computer Science Departmnent, University of Toronto, 1998. |
File: | ge.98.msc.ps.gz |
Title: | Distributed Parallel Shooting for BVODEs |
Authors: | Kwok L. Chow and Wayne H. Enright |
Date: | November 1997 |
Pages: | 57 |
Note: | A revised version of this paper will appear in the ASTC High Performance Computing 1998 |
File: | cs-97-310.ps.gz |
Title: | Fitting Random Stable Solar Systems to Titius-Bode Law |
Authors: | Wayne Hayes and Scott Tremaine |
Date: | October 1997 |
Pages: | 12 |
Note: | Submitted to ICARUS. |
File: | bode.97.ps.gz |
Title: | Superconvergent Interpolants for the Collocation Solution of BVODEs |
Authors: | W. H. Enright and P. H. Muir |
Date: | October 21, 1997 |
Pages: | 42 |
Note: | Address of P. H. Muir: Department of Mathematics and Computing Science Saint Mary's University, Halifax, Nova Scotia, Canada B3H 3C3 Paul.Muir@stmarys.ca |
File: | colintrp97.ps.gz |
Title: | The Parallel Solution of Almost Block Diagonal Systems Arising in Numerical Methods for BVPs for ODEs |
Authors: | R. N. Pancer and K. R. Jackson |
Date: | March 1997 |
Pages: | 6 |
Note: | In the Proceedings of the Fifteenth IMACS World Congress, Berlin, Germany, 1997, Vol. 2, pp. 57-62. |
File: | imacs.97.ps.gz |
Title: | Multigrid and Multilevel Methods for Quadratic Spline Collocation |
Authors: | Christina C. Christara and Barry Smith |
Date: | February 1997 |
Pages: | 21 |
Note: | Appeared in BIT 37:4 (1997), pp. 781-803. |
File: | multig.ps.gz |
DOI: | 10.1007/BF02510352 |
Journal: | PDF file (fulltext), page |
Title: | Validated Solutions of Initial Value Problems for Ordinary Differential Equations |
Authors: | N. S. Nedialkov, K. R. Jackson and G. F. Corliss |
Date: | February 1997 |
Pages: | 40 |
Note: | A slightly revised version of this paper appeared in the journal Applied Mathematics and Computation, 105 (1999), pp. 21-68. |
File: | vsode.97.ps.gz |
Title: | Scientific Computing by Numerical Methods |
Author(s): | C. C. Christara and K. R. Jackson |
Date: | 1996 |
Pages: | 121 |
File: | survey.pdf |
Note: | A revised version of this article appeared in Encyclopaedia of Applied Physics, Vol 17, 1996, pp. 1-79, and in Mathematical Tools for Physicists, 2005, pp. 281-383. |
Title: | Global Error Measures for Large N-body Simulations |
Authors: | Wayne B. Hayes and Kenneth R. Jackson |
Date: | October 1996 |
Pages: | 3 |
Note: | Accepted for publication in the Proceedings of the 12th `Kingston meeting' on Theoretical Astrophysics: Computational Astrophysics, ASP Conference Series, Vol. ??, eds. D. A. Clarke & M. J. West, (San Francisco: ASP). |
File: | shadowing-96b.ps.gz |
Title: | DIMSEMs -- Diagonally Implicit Single-Eigenvalue Methods for the Numerical Solution of Stiff ODEs on Parallel Computers |
Authors: | R. F. Enenkel and K. R. Jackson |
Date: | September 1996 |
Pages: | 32 |
Note: | A slightly revised version of this paper appeared in the journal Advances in Computational Mathematics, 7 (1997), pp. 97-133. |
File: | dimsem.96.ps.gz |
Title: | Remarks on the optimal convolution kernel for CSOR waveform relaxation |
Authors: | Min Hu, Ken Jackson, Jan Janssen and Stefan Vandewalle |
Date: | August 1996 |
Pages: | 24 |
Note: | A slightly revised version of this paper appeared in the journal Advances in Computational Mathematics, 7 (1997), pp. 135-156. |
File: | csor.96.ps.gz |
Title: | Fast Fourier transform solvers for quadratic spline collocation |
Author: | A. Constas |
Date: | July 1996 |
Pages: | 64 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 1996. |
File: | ccc/constas-96-msc.ps.gz |
Title: | DIMSEMs -- Diagonally Implicit Single-Eigenvalue Methods for the Numerical Solution of Stiff ODEs on Parallel Computers |
Author: | Robert F. Enenkel |
Date: | 1996 |
Pages: | 199 |
Note: | Ph.D. Thesis, Computer Science Dept., Univ. of Toronto, 1996. |
File: | enenkel-96-phd.ps.gz |
Title: | Fast Shadowing Algorithm for High Dimensional ODE Systems |
Author: | Wayne Hayes and Kenneth R. Jackson |
Date: | March 1996 |
Pages: | 22 |
File: | shadowing-96a.ps.gz |
Conference: | Computational Finance Workshop |
Organizers: | Peter Forsyth, Computer Science, University of Waterloo Ken Jackson, Computer Science, University of Toronto |
Date: | Thursday, May 2, 1996 |
Note: | Files associated with the Computational Finance Workshop held at the University of Toronto, May 2, 1996 |
Directory: | comp.fin.workshop |
File: | comp.fin.workshop/agenda |
Title: | Modeling Path Dependent Options |
Author: | Yidong Liu |
Date: | May 1996 |
Pages: | 20 |
Note: | Computational Finance Workshop, University of Toronto, May 2, 1996 |
File: | comp.fin.workshop/liu.ps.gz |
Title: | PDE Models of the Term -- Structure: Matching Yield and Volatility Curves |
Author: | Ken Vetzal |
Date: | May 1996 |
Pages: | 17 |
Note: | Computational Finance Workshop, University of Toronto, May 2, 1996 |
File: | comp.fin.workshop/vetzal.ps.gz |
Title: | Asset Allocation with Transaction Costs |
Author: | Hailiang Yang |
Date: | May 1996 |
Pages: | 22 |
Note: | Computational Finance Workshop, University of Toronto, May 2, 1996 |
File: | comp.fin.workshop/yang.ps.gz |
Title: | Robust Numerical Methods for PDE Models of Asian Options |
Author: | Rob Zvan |
Date: | May 1996 |
Pages: | 14 (zvan.1.ps) and 8 (zvan.2.ps) |
Note: | Computational Finance Workshop, University of Toronto, May 2, 1996 |
Files: | comp.fin.workshop/zvan.1.ps.gz and comp.fin.workshop/zvan.2.ps.gz |
Title: | Numerical Pricing of Path-Dependent Options |
Author: | Yidong Liu |
Date: | April 1996 |
Pages: | 68 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 1996. |
File: | liu-96-msc.ps.gz |
In this thesis we first review basic option pricing theory and several popular numerical methods for option pricing. After discussing the analytic Black-Scholes formula for option prices, we focus on numerical methods for option pricing, including Monte Carlo methods, Tree methods and Finite Difference methods.
The main result of this thesis is an extension of the Hull-White tree method for pricing Asian Options. The original Hull-White method is applicable to options in which the averaging period is taken from the beginning to the end of the life of the option. We relax this restriction to any averaging period that ends with the life of the option. The Hull-White tree method is a special case of our method.
Title: | Runge-Kutta Research at Toronto |
Authors: | T. E. Hull, W. H. Enright and K. R. Jackson |
Date: | May 1996 (replaces an earlier draft of March 1996) |
Pages: | 16 |
Note: | A slightly revised version of this paper appeared in the journal Applied Numerical Mathematics, 22 (1996), pp. 225-236. |
File: | rk.survey.96.ps.gz |
There are several main themes. New Runge-Kutta formulas and new error control strategies are developed, leading for example to continuous methods and to new application areas such as delay, differential- algebraic and boundary value problems. Software design and implementation are also emphasized. And so is the importance of careful testing and comparing. Other topics, such as the notion of effectiveness, taking advantage of parallelism, and handling discontinuities are also discussed.
Title: | Improving the Efficiency of Runge-Kutta Methods for the Solution of BVPs for Higher-Order ODEs |
Author: | Khalid Zuberi |
Date: | January 1996 |
Pages: | 96 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 1996. |
File: | zuberi-96-msc.ps.gz |
Title: | Numerical Solution of Retarded and Neutral Delay Differential Equations using Continuous Runge-Kutta Methods |
Author: | Hiroshi Hayashi |
Date: | 1996 |
Pages: | 78 |
Note: | Ph.D Thesis, Computer Science Dept., Univ. of Toronto, 1996. |
File: | hayashi-96-phd.ps.gz |
Title: | Interpolation and Error Control Schemes for Algebraic Differential Equations Using Continuous Implicit Runge-Kutta Methods |
Author: | Hung Nguyen |
Date: | 1995 |
Pages: | 97 |
Note: | Ph.D Thesis, Computer Science Dept., Univ. of Toronto, 1995. |
File: | nguyen-95-phd.ps.gz |
Title: | The numerical solution of large systems of stiff IVPs for ODEs |
Author: | K. R. Jackson |
Date: | September 1995 |
Pages: | 17 |
Note: | A slightly revised version of this paper appeared in the journal Applied Numerical Mathematics, 20 (1996), pp. 5-20. |
File: | mol.95.ps.gz |
Title: | Improving performance when solving high order and mixed order boundary value problems in ODEs |
Authors: | Wayne H. Enright and Min Hu |
Date: | March 1995 |
Pages: | 9 |
File: | high-order-BVP.ps.gz |
Title: | Efficient Shadowing of High Dimensional Chaotic Systems with the Large Astrophysical N-body Problem as an Example |
Author: | Wayne Hayes |
Date: | January 1995 |
Pages: | 61 (thesis.ps) and 2 (errata.ps) |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 1995. |
Directory: | hayes-95-msc |
Files: | hayes-95-msc/thesis.ps.gz and hayes-95-msc/errata.ps.gz |
Title: | DIMSEMs -- Diagonally Implicit Single-Eigenvalue Methods for the Numerical Solution of Stiff ODEs on Parallel Computers: A Preliminary Report |
Author: | R. F. Enenkel and K. R. Jackson |
Date: | January 1995 |
Pages: | 138 |
File: | dimsem.prelim.ps.gz |
Title: | Parallel solvers for spline collocation equations |
Author: | C. C. Christara |
Date: | January 1995 |
Pages: | 35 |
Note: | Appeared in Advances in Engineering Software, 27:1/2, 1996, pp. 71-89 |
File: | ccc/civil.ps.gz |
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: | Precision Control and Exception Handling in Scientific Computing |
Author: | K. R. Jackson and N. S. Nedialkov |
Date: | 1995 |
Pages: | 8 |
File: | prec.except.ps.gz |
The proposed precision control and exception handling facilities are illustrated by sample SciLib programs.
Title: | Precision Control and Exception Handling in Scientific Computing |
Author: | Nedialko (Ned) Stoyanov Nedialkov |
Date: | September 1994 |
Pages: | 51 |
Note: | M.Sc. Thesis, Computer Science Dept., Univ. of Toronto, 1994. |
File: | ned-94-msc.ps.gz |
Precision of computations can be changed during program execution. The first two precisions are IEEE single and double. Increasing the precision provides at least two times more significant digits while the exponent range is significantly larger than double the exponent range of the previous precision.
The exception handling mechanism treats only numerical exceptions and does not distinguish between different types of numerical exceptions.
A variable-precision and exception handling library, SciLib, has been implemented in C++. A new scalar data type real is introduced consisting of variable-precision floating-point numbers. Arithmetic, relational, input and output operators of the language are overloaded for the real data type, so it can be used like any other scalar data type of C++.
The proposed precision control and exception handling are illustrated using SciLib. The examples show the benefits of using precision control to handle exceptions and to avoid catastrophic cancelations.
Title: | Interpolating Runge-Kutta methods for vanishing delay differential equations |
Authors: | Wayne H. Enright and Min Hu |
Date: | 1994 |
Pages: | 34 |
Ref: | Technical report 292/94, Dept. Computer Science, Univ. of Toronto |
Note: | The files int_ddrk6n.f.gz and int_ddvss6.f.gz contain FORTRAN format of the coefficients (in Appendix A and B of the report) of the special interpolants for the codes DDRK6N and DDVSS6 respectively. |
Files: | cs-94-292.ps.gz, int_ddrk6n.f.gz and int_ddvss6.f.gz |
Title: | An Application of Runge-Kutta Predictor-Corrector Methods To Two System of Hyperbolic Equations Arising in Optics |
Author: | K. R. Jackson |
Date: | July 1994 |
Pages: | 3 |
Note: | Appeared in the Proceedings of the Fourteenth IMACS World Congress, Atlanta GA, July, 1994, W. F. Ames ed., pp.~1329--1331 |
File: | imacs.94.ps.gz |
Title: | A Survey of the Explicit Runge-Kutta Method |
Authors: | Wayne H. Enright, Desmond J. Higham, Brynjulf Owren and Philip W. Sharp |
Date: | April 1995 |
Pages: | 36 |
Ref: | Technical report 291/94, Dept. Computer Science, Univ. of Toronto |
Note: | This is a REVISED version, dated April 1995 |
File: | cs-94-291.ps.gz |
We describe recent work in the derivation of Runge-Kutta coefficients: ``classical'' general-purpose formulas, ``special'' formulas for high order and Hamiltonian problems, and ``continuous'' formulas for dense output. We also give a thorough review of implementation details. Modern techniques are described for the tasks of controlling the local error in a step-by-step integration, computing reliable estimates of the global error, detecting stiffness, and detecting and handling discontinuities and singularities. We also discuss some important software issues.
Title: | High performance computing of elliptic partial differential equations with spline collocation |
Author: | C. C. Christara |
Date: | June 1994 |
Pages: | 13 |
Note: | Appeared in Proceedings of the 1994 Computational Structures Technology conference (CST94), August 30 - September 1, 1994, Athens, Greece, sponsored by Civil-Comp, Vol. A, pp. 23-34 |
File: | ccc/highperf.ps.gz |
Title: | An Analysis of the Order of Runge-Kutta Methods That Use an Iterative Scheme To Compute Their Internal Stage Values |
Authors: | K. R. Jackson, A. Kvaerno and S. P. Norsett |
Date: | May 1994 |
Pages: | 53 |
Note: | A slightly revised version of this paper appeared in the journal BIT, 36 (1996), pp. 713-765. |
File: | tr-94-1.ps.gz |
Title: | Quadratic spline collocation methods for elliptic partial differential equations |
Author: | C. C. Christara |
Date: | March 1994 (replaces Tech. Rep. DCS-135, 1990) |
Pages: | 34 |
Note: | A slightly revised version of this paper appeared in BIT, 34:1, 1994, pp. 33-61 |
File: | ccc/cs-90-135.ps.gz |
DOI: | 10.1007/BF01935015 |
Journal: | PDF file (fulltext), page |
Title: | The Mathematical Basis for a New Polynomial Rootfinder with Quadratic Convergence |
Authors: | T. E. Hull and R. Mathon |
Date: | 1994 |
Pages: | 16 |
File: | hull-mathon-math-basis.ps.gz |
This paper first provides the relevant mathematical results:a short derivation of the formulas, convergence proofs, an indication of the behavior near a multiple root, and finally some error bounds. It then provides the outline of an algorithm based on these results, along with some preliminary experimental results to illustrate the major points, and to demonstrate some of the potential for a future software package. It should also be noted that the method is well-suited to take advantage of a parallel environment.
Title: | Coupled Mode Equations with Free Carrier Effects: A Numerical Solution |
Authors: | N. G. R. Broderick, C. M. de Sterke and K. R. Jackson |
Date: | May 1993 |
Pages: | 15 |
Note: | A slightly revised version of this paper appeared in the J. Optical and Quantum Electronics, 26 (1994), pp. S219-S234. |
File: | joqe-1.ps.gz |
Title: | The Use of Butcher Series in the Analysis of Newton-Like Iterations in Runge-Kutta Formulas |
Authors: | K. R. Jackson, A. Kvaerno and S. P. Norsett |
Date: | April 1993 |
Pages: | 17 |
Note: | A slightly revised version of this paper appeared in the journal Applied Numerical Mathematics, 15 (1994), pp. 341-356. |
File: | anm-1.ps.gz |
Title: | A Runge-Kutta Type Boundary Value ODE Solver with Defect Control |
Authors: | W.H. Enright and Paul Muir |
Date: | June 1993 |
Pages: | 41 |
File: | cs-93-267.ps.gz |
Title: | Implementing Complex Elementary Functions Using Exception Handling |
Authors: | T. E. Hull, Thomas F. Fairgrieve and Ping Tak Peter Tang |
Date: | 1993 |
Pages: | 31 |
File: | cmplx-elem-fcns.ps.gz |
Title: | The Practical Implications of Order Reduction on the Solution of Stiff Initial Value Problems using Runge-Kutta Formulas. |
Authors: | C. MacDonald and W. Enright |
Date: | June 1993 |
Pages: | 14 |
Ref: | Technical report 261/92, Dept. Computer Science, Univ. of Toronto |
File: | cs-92-261.ps.gz |
In this paper we examine the phenomenon of `order reduction' and its implications on variable-step algorithms. We introduce a global measure of order referred to here as the `observed order' which is based on the average stepsize over the region of integration. This measure may be better suited to the study of stiff systems, where the stepsize selection algorithm will vary the stepsize considerably over the interval of integration. Observed order gives a better indication of the relationship between accuracy and cost. Using this measure, the `observed order reduction' will be seen to be less severe than that predicted by fixed stepsize order analysis.
Title: | Parallel computation of partial differential equations on distributed memory machines |
Author: | C. C. Christara |
Date: | June 1992 |
Pages: | 13 |
Note: | Appeared in Advances in Computer Methods for Partial Differential Equations - VII, pp. 142-149 R. Vichnevetsky, D. Knight, G. Richter eds. |
File: | ccc/paral.ps.gz |
Title: | The Parallel Solution of ABD Systems Arising in Numerical Methods for BVPs for ODEs |
Authors: | K. R. Jackson and R. N. Pancer |
Date: | May 1992 |
Pages: | 35 |
Ref: | Technical report 255/91, Dept. Computer Science, Univ. of Toronto |
File: | cs-91-255.ps.gz |
In recent papers, Wright presents new stable parallel algorithms for solving ABD systems. Each algorithm attains the theoretically optimal speedup for the problem if enough processors are available, and each can be adapted for use on architectures with fewer than the optimal number of processors. The algorithms in Section 3 of this paper were discovered independently and are similar to those investigated by Wright. Extensions to the basic algorithm are described which make better use of available processors during the decomposition stage in order to increase parallelism in the back-solve stage.
A new algorithm based on "eigenvalue rescaling" is presented in Section 4 which can be up to twice as fast as the fastest algorithm in Section 3. Numerical experiments show that this new algorithm does not exhibit the instability inherent in some earlier schemes. In addition, some insight supporting the numerical evidence is given.
Title: | Order Barriers and Characterizations for Continuous Mono-Implicit Runge-Kutta Schemes |
Authors: | Paul Muir and Brynjulf Owren |
Date: | December 1991 |
Pages: | 22 |
Ref: | Technical report 258/91, Dept. Computer Science, Univ. of Toronto |
File: | cs-91-258.ps.gz |
Title: | Adaptive Linear Equation Solvers in Codes for Large Stiff Systems of ODEs |
Authors: | K. R. Jackson and W. L. Seward |
Date: | July 1991 |
Pages: | 23 |
Ref: | Research Report CS-91-33, Department of Computer Science, University of Waterloo, Waterloo, Ontario, Canada, N2L 3G1 |
Note: | A slightly revised version of this paper appeared in the SIAM J. Sci. Comput., 14 (1993), pp. 800-823. |
File: | uw-cs-91-33.ps.gz |
Title: | Multicolour orderings and iterative methods for elliptic equations |
Author: | C. C. Christara |
Date: | June 1991 |
Pages: | 11 |
Note: | Appeared in Proceedings of the Supercomputing Symposium '91 (SS '91), June 3-5, 1991, Fredericton, NB, Canada, pp. 221-230 |
File: | ccc/colour.ps.gz |
Title: | Domain decomposition and incomplete factorisation methods for partial differential equations |
Author: | C. C. Christara |
Date: | May 1991 |
Pages: | 11 |
Note: | Appeared in Proceedings of the sixth Distributed Memory Computing Conference (DMCC6), April 28 - May 1, 1991, Portland, OR, U.S.A., pp. 369-377 |
File: | ccc/iccg.ps.gz |
Title: | Nonlinear coupled mode equations on a finite interval: a numerical procedure |
Authors: | C. M. de Sterke, K. R. Jackson and B. D. Robert |
Date: | February 1991 |
Pages: | 31 |
Note: | A slightly revised version of this paper appeared in the J. Optical and Quantum Electronics, 26 (1994), pp. S219-S234. |
File: | osa-1.ps.gz |
Title: | A Survey of Parallel Numerical Methods for Initial Value Problems for Ordinary Differential Equations |
Author: | K. R. Jackson |
Date: | 1991 |
Pages: | 6 |
Note: | A slightly revised version of this paper appeared in the IEEE Transactions on Magnetics, 27 (1991), pp. 3792-3797. |
File: | ieee-trans-mag-1.ps.gz |
Title: | Schur complement preconditioned conjugate gradient methods for spline collocation equations |
Author: | C. C. Christara |
Date: | May 1991 |
Pages: | 14 |
Note: | An earlier version of this paper appeared in Computer architecture news, 18:3, 1990, pp. 108-120, and Proceedings of the 1990 International Conference on Supercomputing, June 1990, Amsterdam, the Netherlands, sponsored by ACM |
File: | ccc/scmpcg.ps.gz |
Title: | Conjugate gradient methods for spline collocation equations |
Author: | C. C. Christara |
Date: | May 1991 |
Pages: | 9 |
Note: | An earlier version of this paper appeared in Proceedings of the fifth Distributed Memory Computing Conference (DMCC5), April 8-12, 1990, Charleston, SC, U.S.A., pp. 550-558 |
File: | ccc/pcgcoef.ps.gz |
Title: | Derivation of Efficient Continuous Explicit Runge-Kutta Methods |
Authors: | B. Owren and M. Zennaro |
Date: | 1990 |
Pages: | 17 |
Ref: | Technical report 240/90, Dept. Computer Science, Univ. of Toronto |
File: | cs-90-240.ps.gz |
Title: | The Potential for Parallelism in Runge-Kutta Methods -- Part 1: RK Formulas in Standard Form |
Authors: | K. R. Jackson and S. P. Norsett |
Date: | November 1990 |
Pages: | 36 |
Ref: | Technical report 239/90, Dept. Computer Science, Univ. of Toronto |
Note: | A slightly revised version of this paper appeared in the SIAM J. Numer. Anal., 32 (1995), pp. 49-82. |
File: | cs-90-239.ps.gz |
Title: | Using the IMSL MATH/LIBRARY in Numerical Methods Courses |
Authors: | K. R. Jackson and T. E. Hull |
Date: | April 1990 |
Pages: | 13 |
Note: | A slightly revised version of this paper appeared in the IMSL Technical Report Series, No. 9005. |
File: | imsl-9005.ps.gz |