
SCIENTIFIC PROGRAMS AND ACTIVITIES 

January 31, 2015  
Numerical and Computational Challenges in Science and EngineeringSouthern Ontario Numerical Analysis Day

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
Lowrisk 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 BlackScholes 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 timevarying 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 integralalgebraic 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) functionaldifferential and Volterratype integrodifferential equations with variable delays. Part of the talk has been motivated by certain "integralalgebraic" equations arising in, e.g., memory kernel identification problems in heat conduction or in viscoelasticity, or in the mathematical modeling of lowpressure chemical vapor deposition reactors for growing superconducting films in cellular communications. Since the setting of such IAEs is infinitedimensional (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
ProblemIndependant 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 3D 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 finitedifference 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 problemindependent 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 problemspecific 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
(3D open problems and guidedwave 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 3D 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 twolayer antireflection 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 antireflection 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 104 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.
References:
[1] Yotka Rickard, N. Georgieva and W.P. Huang, “A perfectly matched
layer for the 3D 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 3D wave equation in the time domain”,
in print, IEEE Trans. On Antennas and Propagation, Dec 2002.
[3] N. Georgieva and E. Yamashita, “Timedomain vectorpotential
analysis of transmission line problems,” IEEE Trans. Microwave
Theory Tech., vol. 46, No. 4, 1998, pp. 404410.
[4] N. Georgieva and Y. Rickard, “The application of the wave potential
functions to the analysis of transient electromagnetic fields”,
IEEE MTTS Int. Symposium Digest, June 2000, vol. 2, pp. 11291132.
[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.599604.
[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. 302307
[7] J.P. Berenger, “A perfectly matched layer for the absorption
of electromagnetic waves,” J. of Comp. Physics, vol. 114, 1994,
pp. 185200.
[8] B. Chen, D.G. Fang and B.H. Zhou, “Modified Berenger PML absorbing
boundary condition for FDTD meshes,” IEEE Microwave and Guided
Wave Letters, vol.5, No. 11, 1995, pp. 399401.
[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.
451453.
[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. 914916
[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. 15061509.
[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. 10551063.
[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.25355.
[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 finitedifference timedomain
method,” IEEE Photonics Technology Letters, vol.13, No.5, 2001,
pp. 454456.
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 onesided Lipschitz (OSL) condition.
First, I will consider finite time, strong convergence. I will show
that an implicit variant of the EulerMaruyama 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 longtime asymptotic property:
the ability of numerical methods to reproduce exponential mean square
stability of SDEs. A simple counterexample demonstrates that, for nonglobally
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 longtime 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 NonUniform
Partitions
Quadratic and Cubic Spline Collocation (QSC and CSC) methods of optimal orders of convergence have been developed for the solution of linear secondorder twopoint Boundary Value Problems discretized on uniform partitions. In this paper, we extend the optimal QSC and CSC methods to nonuniform partitions. To do this, we use a mapping function between uniform and nonuniform partitions and develop expansions of the error at the nonuniform 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 (h3j) globally and O (h4j) locally on certain points. The jth derivative of the CSC approximation, is O (h4j) globally, for j > 0, and O (h5j) locally on certain points, for j > 0. The nonuniform 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, ccc@cs.utoronto.ca
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 (nonlinear) 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 GonzalezVega
and Dr. Gema DiazToca of the University of Cantabria, Spain, have been
very helpful.
P.W.Sharp, University of Auckland
Long simulations of the Solar System: symplectic or nonsymplectic?
Long simulations of the Solar System require integrations spanning
millions of years. If the simulations are done using nonsymplectic
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 nonsymplectic methods. Roundoff error complicates
the choice. I will present numerical comparisons between nonsymplectic
and symplectic RungeKutta Nyström methods on several realistic
simulations of the Solar System and discuss the tradeoffs 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 SelfRegular Interior Point Methods
Coauthors: J. Peng, T. Terlaky, G. Zhang
Largeupdate IPMs, induced by Selfregular proximities, form a new class
of primaldual IPMs. These IPMs enjoy better theoretical worst case
complexity than traditional largeupdate IPMs. In this talk we will
discuss our experiences with a selfregular 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, predictorcorrector, dynamic step length, and numerical stability
are addressed and evaluated. The results of our numerical tests on IBM
RS6000 workstations are reported and compared with the LIPSOL implementation,
IBM OSL's implementation, and our implementation that is based on a
selfregular infeasible model.
Attila M. Zsaki, University of Toronto
Computational Challenges in LargeModel 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.