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.
2009
2008
2004
2001
2000
1999
1998
1997
1996
1995
1994
1993
1992
1991
1990
| Title: | On the Blunting Method in the Verified Integration of ODEs |
| Authors: | N. S. Nedialkov, K. R. Jackson and M. Neher |
| Date: | May 2009 |
| Pages: | 17 |
| Note: | Submitted to the Proceedings of the Fifth International Workshop on Taylor Model Methods |
| File: | ned_jac_neh_blunting_2009a.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 and maintaining a low level of overestimation at the same time.
This paper is concerned with one of the sources of overestimation: a
matrix-vector product describing a parallelepiped in the phase space. We
analyze the so-called 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) = u0 ∈ u0,
where u ∈ Rm, A ∈ Rm x m, m ≥ 2
and u0 is a given interval vector, we compare
the convergence behavior of the blunting method with that of the well-known QR
interval method.
For 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: | Randomization in the First Hitting Time Problem |
| Authors: | Ken Jackson, Alex Kreinin and Wanhe Zhang |
| Date: | February 2009 |
| Pages: | 10 |
| Note: |
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. The latest version is available below. |
| 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: | 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 will appear in the International Journal of Theoretical and Applied Finance. |
| 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: | 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: | to appear in ACM Trans. on Math. Soft. |
| File: | enm07.ps.gz |
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: | On Exponential Approximation to the Hockey Stick Function |
| Authors: | Ian Iscoe, Ken Jackson, Alex Kreinin and Xiaofang Ma |
| Date: | January 2007 |
| Pages: | 19 |
| Note: | Submitted to the Journal of Applied Numerical Mathematics |
| File: | IJKM.paper2.pdf |
The hockey stick function is a basic function in pricing and risk management of many financial derivatives. This paper considers approximating the hockey stick function by a sum of exponentials. The algorithm proposed by Beylkin and Monzon is used to determine the parameters of an approximation. Theoretical properties of the approximation are studied. Numerical results are presented.
| Title: | Pricing Correlation-Dependent Derivatives Based on Exponential Approximations to the Hockey Stick |
| Authors: | Ian Iscoe, Ken Jackson, Alex Kreinin and Xiaofang Ma |
| Date: | January 2007 |
| Pages: | 19 |
| Note: | Submitted to the Journal of Computational Finance |
| File: | IJKM.paper3.pdf |
Correlation-dependent derivatives, such as Asset-Backed Securities (ABS) and Collateralized Debt Obligations (CDO), have grown rapidly. 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.ps.gz |
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: | 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: | to appear (2005) in Applied Numerical Mathematics (APNUM) accepted by eds., March 2005. |
| File: | enr05b.ps.gz |
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: | to appear (2005) in Journal of Computational and Applied Mathematics (JCAM) accepted by eds., February 2005. |
| File: | heso05.ps.gz |
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: | to appear in Journal of Nonlinear Analysis, 2005. Also appeared in the Proceedings of the Fourth World Congress of Nonlinear Analysis, Orlando, June 2004. |
| File: | hesg05.ps.gz |
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: | Submitted to the Journal of Risk |
| 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: | 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:
| 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 |