January 21, 2019

Numerical and Computational Challenges in Science and Engineering

Southern Ontario Numerical Analysis Day
Friday, April 19, 2002

Speaker Abstracts

Amélie Bélanger Yotka Rickard
Hermann Brunner Azar Shakoori
Carol Cheng P.W.Sharp
Natalia Georgieva Wen Zhang
D J Higham Lois Zhu
Kit Sun Ng Attila M. Zsaki

Amélie Bélanger, University of Waterloo
Computing the Deltas - Mesh generation for accurate Gradient Computation

Low-risk management of options to buy or sell assets depends on the ability to effectively control the risk associated with such derivatives through hedging. Under ideal Black-Scholes model assumptions, the price of an option and a corresponding hedging strategy (referred to as delta hedging) can be computed from the solution of PDE's. The ideal hedging strategy requires the continual computation of the spatial derivatives of this solution.

We consider the hedging of a European call option written on the max of two underlying assets. Such a problem is studied on a twodimensional domain as the value of the option depends on the value of both underlying assets. Consequently, the hedging strategy in this case will depend on both components of the gradient. We wish to design an effective mesh that will provide acceptable gradient accuracy over the whole domain considered.

A second element that must be taken into account during mesh design is the time-varying nature of the problem. Since the value of the option is a function of time, the solution profile will become sharper as we near maturity date. Our mesh was hence designed to accommodate for this changing solution profile.

Preliminary error analysis with the designed mesh was carried out on a simple surface with a hyperbolic profile intended to represent the payoff at a given time before expiry. Similar error analysis was carried out using numerical data obtained when solving the PDE equation numerically. Finally, simulations will be carried out to confirm that our mesh provides accurate derivative estimates throughout the life of an option.

Hermann Brunner, Memorial University of Newfoundland
From functional differential equations to integral-algebraic equations: theoretical and computational aspects of collocation methods

In this talk I shall describe, in survey style, recent results and a number of open problems in the numerical analysis of collocation methods for (systems of) functional-differential and Volterra-type integro-differential equations with variable delays. Part of the talk has been motivated by certain "integral-algebraic" equations arising in, e.g., memory kernel identification problems in heat conduction or in viscoelasticity, or in the mathematical modeling of low-pressure chemical vapor deposition reactors for growing super-conducting films in cellular communications. Since the setting of such IAEs is infinite-dimensional (they may be viewed as abstract DAEs), new approaches for their numerical analysis and computational solution have to be explored. Thus, the current research in this area I am involved in with colleagues at Humboldt University and Virginia Tech encompasses a wide spectrum of mathematics and scientific computation and contains an equally wide range of challenging open problems.

Carol Cheng, University of Waterloo
Demand Deposit Valuation

Joint work with Peter Forsyth, University of Waterloo and Ken Vestal, University of Waterloo.

In recent years, core deposits such as checking accounts typically comprise about 35% of commercial bank liabilities. Valuing these liabilities has a large impact on bank balance sheets. This paper builds on previous work on demand deposit valuation. We value demand deposits based on data from RBC Finance Group, and determine the best numerical scheme under different scenarios. Numerical schemes such as Partial Differential Equation (PDE), Monte Carlo methods including crude Monte Carlo method and Antithetic Monte Carlo method, and Low discrepancy sequence methods (LDS) are tested. Convergence of different numerical schemes is compared. It is shown that the Monte Carlo method is adequate for short time horizons with low accuracy requirements, while the PDE method generates better results in other cases. Our result shows that these deposits are extremely profitable.

Keywords: Demand deposit, PDE, Monte Carlo, LDS

Natalia Georgieva and Yotka Rickard, McMaster University
Problem-Independant Enhancement of PML ABC For Finite Difference Time Domain Techniques in Electrodynamics

Recently, an enhanced perfectly matched layer (PML) absorbing boundary condition (ABC) has been proposed [1], [2] for the 3-D wave equation in the time domain (WETD) derived for the vector potentials (VPs) [3], [4]. The derivation is based on the stretched coordinate approach [5]. Thus the PML is defined in the frequency domain and then six auxiliary variables are introduced to facilitate its mapping into the time domain. In [2], a new degree of freedom in the definition of the PML variable profiles has been introduced in the existing spatially polynomial scaling, and it is combined with modifications in the PML terminating walls. That improves the PML performance in applications involving the wave equation. However, an enhanced performance of this PML ABC is also observed when applied with the classical Yee’s finite-difference timedomain (FDTD) method to the solution of Maxwell’s equations [6].

In this work, the performance of this enhanced PML (EPML) ABC, when applied in Yee’s FDTD algorithm, is compared with existing PML ABCs like the original Berenger’s PML [7], the modified PML (MPML) [8] and the generalized PML (GPML) [9]. Also, a comparison of the reflection levels from the EPML ABC is made when one problem is modeled by Yee’s FDTD method and by the WETD method. Finally, the influences of the new degree of freedom and of the modified termination walls on the PML performance is studied and the best choices are singled out both in the case of Yee’s FDTD model and the WETD model.

In [1], [2] it has been proposed firstly the PML conductivity (responsible for the propagating waves attenuation) and the PML loss factor (responsible for the evanescent waves attenuation) to increase at different exponent rates within the PML medium. This in effect adds a new degree of freedom in the definition of the PML variable profiles. In addition, further performance enhancement is achieved by the use of modified PML terminations. Such a problem-independent approach to improve the PML performance is in contrast with the commonly used approach to optimize a single PML variable profile like in [10], [11], [12], and frequently it is done in a problem-specific environment [13], [14].

To illustrate the better performance of the proposed EPML ABC for Yee’s FDTD applications, the following examples are considered: an infinitesimal dipole radiating in open space and a microstrip line. It is shown that the proposed EPML has a superior performance in all major types of problems (3-D open problems and guided-wave problems) in comparison with commonly used PML ABCs such as MPML and GPML: it offers lower reflection levels in a broader frequency band.

Secondly, a comparison of the spectrum of the reflections is done when one and the same 3-D structure is modeled by Yee’s FDTD method and by the WETD method for the magnetic VP. Thus the equivalency of both methods is validated as well as the excellent performance of EPML in both models. The structure is chosen to be an optical buried waveguide terminated by a two-layer anti-reflection coating. The reason for this choice is that in such a problem the low reflection level from the ABC itself is of utmost importance since the reflections from the anti-reflection coating are of the order of the reflections from most of the commonly used ABCs. By initially terminating an infinite (longitudinally uniform) optical buried waveguide directly by EPML, it is shown that the level of the reflections caused by the EPML itself in both models is below –80 dB (corresponding to 10-4 magnitude).

Finally, the influence of the difference between the exponent rates of the PML conductivity and the PML loss factor on the overall EPML performance is considered both in the case of Yee’s FDTD model and the WETD model. It is shown that if the PML conductivity increases faster than the PML loss factor, the reflection levels can be decreased by half to full order of magnitude.

[1] Yotka Rickard, N. Georgieva and W.-P. Huang, “A perfectly matched layer for the 3-D wave equation in the time domain”, in print, IEEE Microwave and Wireless Components Letters, May 2002.
[2] Yotka Rickard, N. Georgieva and W.-P. Huang, “Application and optimization of PML ABC for the 3-D wave equation in the time domain”, in print, IEEE Trans. On Antennas and Propagation, Dec 2002.
[3] N. Georgieva and E. Yamashita, “Time-domain vector-potential analysis of transmission line problems,” IEEE Trans. Microwave Theory Tech., vol. 46, No. 4, 1998, pp. 404-410.
[4] N. Georgieva and Y. Rickard, “The application of the wave potential functions to the analysis of transient electromagnetic fields”, IEEE MTT-S Int. Symposium Digest, June 2000, vol. 2, pp. 1129-1132.
[5] W.C. Chew and W.H. Weedon, “A 3D perfectly matched medium from modified Maxwell’s equations with stretched coordinates,” Microwave and Optical Technology Letters, vol.7, No. 13, 1994, pp.599-604.
[6] K.S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media”, IEEE Trans. on Antennas and Propagation, vol. 14, 1966, pp. 302-307
[7] J.P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. of Comp. Physics, vol. 114, 1994, pp. 185-200.
[8] B. Chen, D.G. Fang and B.H. Zhou, “Modified Berenger PML absorbing boundary condition for FD-TD meshes,” IEEE Microwave and Guided Wave Letters, vol.5, No. 11, 1995, pp. 399-401.
[9] J. Fang, Zh. Wu, “Generalized Perfectly Matched Layer – An Extension of Berenger’s Perfectly Matched Layer Boundary Condition”, IEEE Microwave and Guided Wave Letters, vol. 5, No. 12, Dec 1995, pp. 451-453.
[10] G. Lazzi, O.P. Gandhi, “On the Optimal Design of the PML Absorbing Boundary Condition for the FDTD Code”, IEEE Trans. on Antennas and Propagation, vol. 45, no. 5, May 1997, pp. 914-916
[11] E.A. Marengo, C.M. Rappaport and E. L. Miller, “Optimum PML ABC Conductivity Profile in FDTD”, IEEE Trans. on Magnetics, vol. 35, no. 3, May 1999, pp. 1506-1509.
[12] S.C. Winton and C.M. Rappaport, “Specifying PML Conductivities by Considering Numerical Reflection Dependencies”, IEEE Trans. on Antennas and Propagation, vol. 48, no. 7, Jul 2000, pp. 1055-1063.
[13] D. Johnson, C. Furse and A. Tripp, - Application and optimization of the Perfectly Matched Layer Boundary Condition for geophysical simulations, Microwave and Optical Technology Letters, vol. 25, No.4, May 20 2000, pp.253-55.
[14] G. Lazzi, C.M. Furse, O.P. Gandhi, “Optimization and design of conductivity profiles for the PML boundary condition and its application to bioelectromagnetic problems”, IEEE Antennas and Propagation Society International Symposium, Digest, Vol. 1, 1997, pp. 486 –489.
[15] D. Zhou, W.P. Huang, C.L. Xu and D.G. Fang, “The perfectly matched layer boundary condition for scalar finite-difference time-domain method,” IEEE Photonics Technology Letters, vol.13, No.5, 2001, pp. 454-456.

D J Higham, University of Strathclyde, Scotland (Visiting the Fields Institute)
Convergence and Stability of Numerical Methods for Nonlinear Stochastic ODEs

Without requiring any background knowledge of stochastic differential equations (SDEs), I will aim to give an overview of some recent results concerning their numerical simulation. The results were proved in collaboration with Xuerong Mao (Strathclyde) and Andrew Stuart (Warwick), and their purpose is to get around the somewhat restrictive global Lipschitz assumption that appears in the traditional analysis of SDE methods. The key idea, borrowed from the deterministic literature, is to weaken the global Lipschitz condition on the drift coefficient to a one-sided Lipschitz (OSL) condition.

First, I will consider finite time, strong convergence. I will show that an implicit variant of the Euler-Maruyama method converges under a OSL condition on the drift. Further, the optimal rate of convergence can be recovered if the drift coefficient is also assumed to behave like a polynomial. I will then look at a long-time asymptotic property: the ability of numerical methods to reproduce exponential mean square stability of SDEs. A simple counterexample demonstrates that, for non-globally Lipschitz SDEs, this form of stability is not inherited, in general, by numerical methods. However, I will show that under a suitable OSL
condition, positive results may be obtained for two implicit methods.

These results emphasize that for long-time simulation on nonlinear SDEs, the choice of numerical method can be crucial.

Kit Sun Ng, University of Toronto
Optimal Quadratic and Cubic Spline Collocation on Non-Uniform Partitions

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 discretized on uniform partitions. In this paper, we extend the optimal QSC and CSC methods to non-uniform partitions. To do this, we use a mapping function between uniform and non-uniform partitions and develop expansions of the error at the non-uniform collocation points of some appropriately defined spline interpolants. We use the Green's function approach to analyze the proposed methods. The existence and uniqueness of the QSC and CSC approximations are shown, under some conditions. Optimal global and local orders of convergence of the spline approximations and derivatives are derived, similar to those of the respective methods for uniform partitions. The jth derivative of the QSC approximation, for j > 0, is O (h3-j) globally and O (h4-j) locally on certain points. The jth derivative of the CSC approximation, is O (h4-j) globally, for j > 0, and O (h5-j) locally on certain points, for j > 0. The non-uniform partition QSC and CSC methods are integrated with adaptive grid techniques, and applied to the solution of a variety of problems, including problems with interior or boundary layers. The numerical results verify the theoretically expected behaviour of the methods.

Joint work with Christina C. Christara,

Azar Shakoori, The University of Western Ontario
Exploration of Bezout Matrices in Maple for solving Bivariate Polynomials

Systems of polynomial equations occur frequently in applications, especially engineering applications. It is therefore important to have good methods to solve these equations. Three important requirements for a method to be considered useful are, first, that it finds all the solutions: second, that it finds them all accurately; and, finally, that it finds them all quickly. It has recently been observed that, contrary to how students are first taught in linear algebra, it is useful to convert a polynomial problem to an eigenvalue problem. This is true in part because polynomials are sensitive to small changes in their input data, whereas matrices are usually much less sensitive, but conversion to eigenvalue problems is also useful in the multivariate case because good numerical methods are available to solve numerically the generalized eigenproblems that result.
In this talk, I will give a short introduction to the Bezout Matrix and the companion matrix pencil. Finally, I will introduce an efficient method to solve a (non-linear) bivariate polynomial system, based on interpretation of resultants as eigenvalue problems. We show some examples at the end.

This work is under the supervision of Professor Robert M. Corless (Applied Math, University of Western Ontario). Professor Laureano Gonzalez-Vega and Dr. Gema Diaz-Toca of the University of Cantabria, Spain, have been very helpful.

P.W.Sharp, University of Auckland
Long simulations of the Solar System: symplectic or non-symplectic?

Long simulations of the Solar System require integrations spanning millions of years. If the simulations are done using non-symplectic methods, the global error typically grows as t2 and the numerical solution at the end point will have few digits of accuracy. The growth of the global error can be reduced to O (t) by using symplectic methods. This reduction comes at some cost and it is not immediately obvious symplectic methods are better than non-symplectic methods. Round-off error complicates the choice. I will present numerical comparisons between non-symplectic and symplectic Runge-Kutta Nyström methods on several realistic simulations of the Solar System and discuss the trade-offs between the two types of methods.
This was joint work with R. Vaillancourt of the University of Ottawa.

Wen Zhang, Oakland University
Morphological Evolution of a 3D Array of Particles under Surface Diffusion

Joint work with Ian Gladwell, Southern Methodist University

During sintering of particles, surface diffusion forces a movement of particle surface that is governed by the surface Laplacian of curvature. We formulate the surface diffusion equation and boundary conditions in a global coordinate system for tracking surface evolution in a closed packed 3D array of particles. Using the method of lines, we develop a computational strategy and discuss its results for a particle surface movement in the 3D array.

Lois Zhu, McMaster University
Computational Experiences with Self-Regular Interior Point Methods

Coauthors: J. Peng, T. Terlaky, G. Zhang
Large-update IPMs, induced by Self-regular proximities, form a new class of primal-dual IPMs. These IPMs enjoy better theoretical worst case complexity than traditional large-update IPMs. In this talk we will discuss our experiences with a self-regular proximity based IPMs that is implemented on LIPSOL and IBM OSL and using WSMP as a large sparse linear system solver. Various implementation issues, such as initial solution, predictor-corrector, dynamic step length, and numerical stability are addressed and evaluated. The results of our numerical tests on IBM RS-6000 workstations are reported and compared with the LIPSOL implementation, IBM OSL's implementation, and our implementation that is based on a self-regular infeasible model.

Attila M. Zsaki, University of Toronto
Computational Challenges in Large-Model Numerical Analysis of Underground Excavations

The finite and boundary elements methods are often employed in the analysis of underground excavations in mining and geotechnical engineering. More than often the sheer size of the discretized geometry coupled with the material constitutive models pushes the limits of current hardware. Moreover, for the case of the large models the "time to solution" can be measured in days. As part of my research I have formulated a framework of "approximate analysis" of the problems using a simplified geometry to obtain solutions that are representative of the general solution but only accurate in a limited region where the analysis was requested to be accurate. The analysis framework was first developed for 2D scenarios, but now it is expanded to address full 3D domains. In addition to geometry simplification, parallel processing is employed to reduce the computation time.

Return to Workshop Page