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.
2008
2004
2001
2000
1999
1998
1997
1996
1995
1994
1993
1992
1991
1990
| Title: | Randomization in the Default Boundary Problem |
| Authors: | Ken Jackson, Alex Kreinin and Wanhe Zhang |
| Date: | February 2008 |
| Pages: | 9 |
| Note: | Also available from http://www.defaultrisk.com/ |
| File: | First.Hitting.Time.2008.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)= μ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: | Submitted to the International Journal of Theoretical and Applied Finance, December 2007. Also available from http://www.defaultrisk.com/ |
| 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 b