Dynamics of Complex Singularities of Nonlinear PDEs: Analysis and Computation
Accepted to appear in SEMA-SIMAI (special volume devoted to invited talks of ICIAM2019) online
Solutions to nonlinear evolution equations exhibit a wide range of interesting phenomena such as shocks, solitons, recurrence, and blow-up. As an aid to understanding some of these features, the solutions can be viewed as analytic functions of a complex space variable. The dynamics of poles and branch point singularities in the complex plane can often be associated with the aforementioned features of the solution. Some of the computational and analytical results in this area are surveyed here. This includes a first attempt at computing the poles in the famous Zabusky–Kruskal experiment that lead to the discovery of the soliton.
Exponential node clustering at singularities for rational approximation, quadrature, and PDEs
(with Lloyd N Trefethen and Yuji Nakatsukasa)
Numeriesche Mathematik 2021 online
Rational approximations of functions with singularities can converge at a rootexponential rate if the poles are exponentially clustered. We begin by reviewing this effect in minimax, least-squares, and AAA approximations on intervals and complex domains, conformal mapping, and the numerical solution of Laplace, Helmholtz, and biharmonic equations by the “lightning” method. Extensive and wide-ranging numerical experiments are involved. We then present further experiments giving evidence that in all of these applications, it is advantageous to use exponential clustering whose density on a logarithmic scale is not uniform but tapers off linearly to zero near the singularity. We propose a theoretical model of the tapering effect based on the Hermite contour integral and potential theory, which suggests that tapering doubles the rate of convergence. Finally we show that related mathematics applies to the relationship between exponential (not tapered) and doubly exponential (tapered) quadrature formulas. Here it is the Gauss–Takahasi–Mori contour integral that comes into play.
Quadratic Padé Approximation: Numerical Aspects and Applications (with Marco Fasondini, Nick Hale, and René Spoerer)
Computer Research and Modeling, Vol. 11, No. 6, pp. 1017-1031 (2019) online
Padé approximation is a useful tool for extracting singularity information from a power series. A linear Padé approximant is a rational function and can provide estimates of pole and zero locations in the complex plane. A quadratic Padé approximant has square root singularities and can, therefore, provide additional information such as estimates of branch point locations. In this paper, we discuss numerical aspects of computing quadratic Padé approximants as well as some applications. Two algorithms for computing the coefficients in the approximant are discussed: a direct method involving the solution of a linear system (well-known in the mathematics community) and a recursive method (well-known in the physics community). We compare the accuracy of these two methods when implemented in floating-point arithmetic and discuss their pros and cons. In addition, we extend Luke’s perturbation analysis of linear Padé approximation to the quadratic case and identify the problem of spurious branch points in the quadratic approximant, which can cause a significant loss of accuracy. A possible remedy for this problem is suggested by noting that these troublesome points can be identified by the recursive method mentioned above. Another complication with the quadratic approximant arises in choosing the appropriate branch. One possibility, which is to base this choice on the linear approximant, is discussed in connection with an example due to Stahl. It is also known that the quadratic method is capable of providing reasonable approximations on secondary sheets of the Riemann surface, a fact we illustrate here by means of an example. Two concluding applications show the superiority of the quadratic approximant over its linear counterpart: one involving a special function (the Lambert W-function) and the other a nonlinear PDE (the continuation of a solution of the inviscid Burgers equation into the complex plane).
Gauss-Hermite Quadrature for the Bromwich Integral
SIAM Journal of Numerical Analysis, Vol. 57, pp. 2200-2216 (2019) online
Several popular methods for the numerical inversion of the Laplace transform are based on quadrature approximation of the Bromwich integral. A key ingredient is contour deformation, and a variety of contours have been proposed. The choice of quadrature rule, however, has always been the rapidly convergent trapezoidal rule. In this paper the classical Gauss--Hermite rule is proposed as an alternative to the trapezoidal rule. By means of a numerical saddle point approach, contour parameters are computed for the case of a parabolic contour where the singularities of the transform are assumed to be located on the nonpositive real axis. In various numerical experiments it is shown that the Gauss-Hermite method with parameters proposed here converges more rapidly than the best methods based on the trapezoidal rule.
A computational exploration of the McCoy-Tracy-Wu solutions of the third Painlevé equation (with Marco Fasondini and Bengt Fornberg)
Physica D, Vol. 363, pp. 18-43 (2018) online
The method recently developed by the authors for the computation of the multivalued Painlevé transcendents on their Riemann surfaces (Fasondini et al., 2017) is used to explore families of solutions to the third Painlevé equation that were identified by McCoy et al. (1977) and which contain a pole-free sector. Limiting cases, in which the solutions are singular functions of the parameters, are also investigated and it is shown that a particular set of limiting solutions is expressible in terms of special functions. Solutions that are single-valued, logarithmically (infinitely) branched and algebraically branched, with any number of distinct sheets, are encountered. The algebraically branched solutions have multiple pole-free sectors on their Riemann surfaces that are accounted for by using asymptotic formulae and Backlund transformations.
Methods for the computation of the multivalued Painlevé transcendents on their Riemann surfaces (with Marco Fasondini and Bengt Fornberg)
Journal of Computational Physics, Vol. 344, pp. 36-50 (2017) online
We extend the numerical pole field solver (B. Fornberg and J.A.C. Weideman, J. Comput. Phys. 230:5957-5973, 2011) to enable the computation of the multivalued Painlevé transcendents, which are the solutions to the third, fifth and sixth Painlevé equations, on multiple sheets of their Riemann surfaces. We display, for the first time, solutions to these equations on multiple Riemann sheets. We also provide numerical evidence for the existence of solutions to the sixth Painlevé equation that have pole-free sectors, known as tronquée solutions.
Contour Integral Solution of Elliptic PDEs in Cylindrical Domains (with Nicholas Hale)
SIAM Journal of Scientific Computing, Vol. 37, pp. A2630-A2655 (2015) online
The solutions of certain elliptic PDEs can be expressed as contour integrals of Dunford type. In this paper efficient contours and quadrature rules for the approximation of such integrals are proposed. The trapezoidal and midpoint rules are used in combination with a conformal mapping that fully exploits the analyticity of the integrand, leading to rapidly converging quadrature formulas of double exponential type. In addition to optimizing the step size of the quadrature formula, implementation aspects such as the solution of the resulting shifted linear systems at each quadrature node are discussed. Numerical examples involving Laplace's equation in a rectangle, a box, and an annular cylinder are presented. Timing and accuracy comparisons of the various implementations are given.
A Computational Overview of the Solution Space of the Imaginary PII Equation (with Bengt Fornberg)
Physica D, Vol. 309, pp. 108-118 (2015) online
The six Painlevé equations were first formulated about a century ago. Since the 1970's, it has become increasingly recognized that they play a fundamental role in a wide range of physical applications. A recently developed numerical pole field solver (J. Comput. Phys. 230 (2011), 5957-5973) now allows their complete solutions spaces to be surveyed across the complex plane. Following such surveys of the PI , PII and PIV equations, we consider here the case of the imaginary PII equation (the standard PII equation, with a change of sign for its nonlinear term). Solutions to this equation share many features with other classes of Painlevé transcendents, including a rich variety of pole field configurations, with connection formulas linking asymptotic behaviors in different directions of the complex plane.
An Improved Talbot Method for Numerical Laplace Transform Inversion (with Benedict Dingfelder)
Numerical Algorithms, Vol. 68, pp. 167--183 (2015) online
The classical Talbot method for the computation of the inverse Laplace transform is improved for the case where the transform is analytic in the complex plane except for the negative real axis. First, by using a truncated Talbot contour rather than the classical contour that goes to infinity in the left half-plane, faster convergence is achieved. Second, a control mechanism for improving numerical stability is introduced.
Optimal Domain Splitting for Interpolation by Chebyshev Polynomials (with Tobin A Driscoll)
SIAM Journal on Numerical Analysis, Vol. 52, pp. 1913-1927 (2014) online
Polynomial interpolants defined using Chebyshev extreme points as nodes converge uniformly at a geometric rate when sampling a function that is analytic on an interval. However, the convergence rate can be arbitrarily close to unity if the function has a singularity close to the interval when extended to the complex plane. In such cases splitting the interval and doing piecewise interpolation may be more efficient in the total number of nodes than the global interpolant. Because the convergence rate is determined by Bernstein ellipses obtained through a Joukowski conformal map, relative efficiency of splitting at any point in the interval can be calculated and then optimized over the interval. The optimal splitting may be applied recursively. The Chebfun software project uses a simple rule of thumb without prior singularity information to create a binary search that can be shown to do an excellent job of finding the optimal splitting in most cases. However, the process can use a large number of intermediate function evaluations that are not needed in the final approximant. A Chebyshev-Padé technique generates approximate singularity locations that are good enough to get close to the optimal splitting more directly in test cases. The technique is applied to a singularly perturbed boundary-value problem.
The Exponentially Convergent Trapezoidal Rule (with Lloyd N Trefethen)
SIAM Review, Vol. 56, pp. 385-458 (2014) online. Comments and errata.
It is well known that the trapezoidal rule converges geometrically when applied to analytic functions on periodic intervals or the real line. The mathematics and history of this phenomenon are reviewed and it is shown that far from being a curiosity, it is linked with computational methods all across scientific computing, including algorithms related to inverse Laplace transforms, special functions, complex analysis, rational approximation, integral equations, and the computation of functions and eigenvalues of matrices and operators.
A Computational Exploration of the Second Painlevé Equation (with Bengt Fornberg)
Foundations of Computational Mathematics, Vol. 14, pp. 985-1016 (2014)
online
The pole field solver developed recently by the authors (J. Comput. Phys. 230:5957-5973, 2011) is used to survey the space of solutions of the second Painlevé equation that are real on the real axis. This includes well-known solutions such as the Hastings-McLeod and Ablowitz-Segur type of solutions, as well as some novel solutions. The speed and robustness of this pole field solver enable the exploration of pole dynamics in the complex plane as the parameter and initial conditions of the differential equation are varied. Theoretical connection formulas are also verified numerically.
A Numerical Methodology for the Painlevé Equations (with Bengt Fornberg)
Journal of Computational Physics, Vol. 230, pp. 5957-5973 (2011)
The six Painlevé transcendents PI-PVI have both applications and analytic properties that make them stand out from most other classes of special functions. Although they have been the subject of extensive theoretical investigations for about a century, they still have a reputation for being numerically challenging. In particular, their extensive pole fields in the complex plane have often been perceived as "numerical mine fields". In the present work, we note that the Painlevé property in fact provides the opportunity for very fast and accurate numerical solutions throughout such fields. When combining a Taylor/Padé-based ODE initial value solver for the pole fields with a boundary value solver for smooth regions, numerical solutions become available across the full complex plane. We focus here on the numerical methodology, and illustrate it for the PI equation. In later studies, we will concentrate on mathematical aspects of both the PI and the higher Painlevé transcendents.
A Contour Integral Method for the Black-Scholes and Heston Equations
(with Karel in't Hout)
SIAM Journal on Scientific Computing, Vol. 33, pp.763-785 (2011)
A contour integral method recently proposed by Weideman [IMA J. Numer. Anal., 30:334--350, 2010] for integrating semi-discrete advection-diffusion PDEs, is improved and extended for application to some of the important equations of mathematical finance. Using estimates for the numerical range of the spatial operator, optimal contour parameters are derived theoretically and tested numerically. An improvement on the existing method is the use of Krylov methods for the shifted linear systems, the solution of which represent the major computational cost of the algorithm. Test examples presented are the Black--Scholes PDE in one space dimension and the Heston PDE in two dimensions. In the latter case efficiency is compared to ADI splitting schemes for solving this problem. In the experiments it is found that the contour integral method is superior for the range of medium to high accuracy requirements.
Improved Contour Integral Methods for Parabolic PDEs
IMA Journal of Numerical Analysis, Vol. 30, pp. 334-350 (2010)
One way of computing the matrix exponential that arises in semi-discrete parabolic PDEs is via the Dunford-Cauchy integral formula. The integral is approximated by the trapezoidal or midpoint rules on a Hankel contour defined by a suitable change of variables. In a recent paper by the author and Trefethen [Math. Comp., Vol. 76, 1341--1356 (2007)] two widely used contours were analyzed. Estimates for the optimal parameters that define these contours were proposed. In this paper this analysis is extended in two directions. First, the effect of roundoff error is now included in the error model. Second, we extend the results to the case of a model convection-diffusion equation, where a large convective term causes the matrix to be highly non-normal.
High Accuracy Representation of the Free Propagator (with P Nash)
Applied Numerical Mathematics, Vol. 59, pp. 2937--2949 (2009)
A new matrix approximation to the free propagator of quantum mechanics is proposed. It is a variant of the well-known Fourier pseudospectral approximation. Numerical comparisons involving the linear and nonlinear Schr\"odinger equations are presented.
Algorithm 882: Near Best Fixed Pole Rational Interpolation
With Application in Spectral Methods
(with J van Deun,
K Deckers, A Bultheel)
ACM Transactions on Mathematical Software, Vol. 35, pp. 14:1-14:21 (2008)
We present a numerical procedure to compute the nodes and weights in rational Gauss-Chebyshev quadrature formulas. Under certain conditions on the poles, these nodes are near best for rational interpolation with prescribed poles (in the same sense that Chebyshev points are near best for polynomial interpolation). As an illustration, we use these interpolation points to solve a differential equation with an interior boundary layer using a rational spectral method. The algorithm to compute the interpolation points (and, if required, the quadrature weights) is implemented as a MATLAB program.
The Kink Phenomenon in Fejér and Clenshaw-Curtis Quadrature (with LN Trefethen)
Numerische Mathematik, Vol. 107, pp. 707-727 (2007) (online)
The Fejér and Clenshaw-Curtis rules for numerical integration exhibit a curious phenomenon when applied to certain analytic functions. When N (the number of points in the integration rule) increases, the error does not decay to zero evenly but does so in two distinct stages. For N less than a critical value, the error behaves like O(1/r^{2N}), where r is a constant greater than 1. For these values of N the accuracy of both the Fejer and Clenshaw-Curtis rules is almost indistinguishable from that of the more celebrated Gauss-Legendre quadrature rule. For larger N, however, the error decreases at the rate O(1/r^{N}), i.e., only half as fast as before. Convergence curves typically display a kink where the convergence rate cuts in half. In this paper we derive explicit as well as asymptotic error formulas that provide a complete description of this phenomenon.
Parabolic and Hyperbolic Contours for Computing the Bromwich Integral (with LN Trefethen)
Mathematics of Computation, Vol. 76, pp. 1341-1356 (2007) (online)
Some of the most effective methods for the numerical inversion of the Laplace transform are based on the approximation of the Bromwich contour integral. The accuracy of these methods often hinges on a good choice of contour, and several such contours have been proposed in the literature. Here we analyze two recently proposed contours, namely a parabola and a hyperbola. Using a representative model problem, we determine estimates for the optimal parameters that define these contours. An application to a fractional diffusion equation is presented.
Optimizing Talbot's Contours for the Inversion of the Laplace Transform
SIAM Journal on Numerical Analysis, Vol. 44, pp. 2342-2362 (2006) (online)
Talbot's method for the numerical inversion of the Laplace Transform consists of numerically integrating the Bromwich integral on a special contour by means of the trapezoidal or midpoint rules. In this paper we address the issue of parameter selection in the method, for the particular situation when parabolic PDEs are solved. In the process the well-known subgeometric convergence rate O(exp(-c sqrt(N)) of this method is improved to the geometric rate $O(exp(-c N))$, with N the number of nodes in the integration rule. The value of the maximum decay rate c is explicitly determined. Numerical results for two versions of the heat equation are presented. With the choice of parameters derived here, the rule-of-thumb is that to achieve an accuracy of 10^(-ell) at any given time, the associated elliptic problem has to be solved no more than ell times.
Spectral Differentiation Matrices for the Numerical Solution of Schrödinger's Equation
Journal of Physics A: Mathematics and General, Vol. 39, pp. 10229-10237 (2006) (online)
Schrödinger's equation is solved numerically using the spectral collocation (pseudospectral) method based on Hermite weighted polynomial interpolants. Several sets of numerical results from the literature are reproduced using this approach. MATLAB codes that can serve as templates for further exploration are provided both in the paper and online. A new proposal is made for solving Schrödinger's equation in Stokes wedges.
MATLAB materials related to this paper can be found here.
Talbot Quadratures and Rational Approximations (with LN Trefethen and T Schmelzer)
BIT Numerical Mathematics, Vol. 46, pp. 653-670 (2006) (online)
Many computational problems can be solved with the aid of contour integrals containing exp(z) in the the integrand: examples include inverse Laplace transforms, special functions, functions of matrices and operators, parabolic PDEs, and reaction-diffusion equations. One approach to the numerical quadrature of such integrals is to apply the trapezoid rule on a Hankel contour defined by a suitable change of variables. Optimal parameters for three classes of such contours have recently been derived: (a) parabolas, (b) hyperbolas, and (c) cotangent contours, following Talbot in 1979. The convergence rates for these optimized quadrature formulas are very fast: roughly O(3^{-N}), where N is the number of sample points or function evaluations. On the other hand, convergence at a rate apparently about twice as fast, O(9.28903^{-N}), can be achieved by using a different approach: best supremum-norm rational approximants to exp(z) for z in (-infinity,0], following Cody, Meinardus and Varga in 1969. (All these rates are doubled in the case of self-adjoint operators or real integrands.) It is shown that the quadrature formulas can be interpreted as rational approximations and the rational approximations as quadrature formulas, and the strengths and weaknesses of the different approaches are discussed in the light of these connections. A MATLAB function is provided for computing Cody-Meinardus-Varga approximants by the method of Caratheodory-Fejér approximation.
Padé Approximations to the Logarithm I: Derivation via Differential Equations
Quaestiones Mathematicae Vol. 28, pp. 375--390 (2005)
Explicit formulas for the polynomials that define Hermite-Padé approximations to the logarithm are derived. The main elements of the derivation are Frobenius' method for solving linear differential equations, and automatic summation algorithms. Applications to finite difference formulas and Laplace transforms are indicated.
Padé Approximations to the Logarithm II: Identities, Recurrences, and Symbolic Computation (with K Driver, H Prodinger, C Schneider)
Ramanujan Journal, Vol. 11, pp. 139-158 (2006) (online)
Combinatorial identities that were needed in [I] are proved, mostly with C. Schneider's computer algebra package SigmaP. The form of the Padé approximation of the logarithm of arbitrary order is stated as a conjecture.
Padé Approximations to the Logarithm III: Alternative Methods and Additional Results (with K Driver, H Prodinger, C Schneider)
Ramanujan Journal, Vol. 12, pp. 299-314 (2006) (online)
We use C. Schneider's summation software SigmaP and a method due to Andrews-Newton-Zeilberger to reprove results from [I] and [II] as well as to prove new related material.
The Accuracy of the Chebyshev Differencing Method for Analytic Functions (with SC Reddy)
SIAM Journal on Numerical Analysis, Vol. 42, pp. 2176--2187 (2005) (online)
The Chebyshev spectral collocation method is one of the most powerful tools for numerical differentiation, particularly when the function under consideration is smooth. An upper bound on the error, in the discrete maximum norm, is derived here when the function is analytic. Two model functions are analyzed in detail.
Computing the Dynamics of Complex Singularities of Nonlinear PDEs
SIAM Journal on Applied Dynamical Systems, Vol. 2, pp. 171--186 (2003) (electronic)
A two-step strategy is proposed for the computation of singularities in nonlinear PDEs. The first step is the numerical solution of the PDE using a Fourier spectral method; the second step involves numerical analytical continuation into the complex plane using the epsilon algorithm to sum the Fourier series. Test examples include the inviscid Burgers and nonlinear heat equations, as well as a transport equation involving the Hilbert transform. Numerical results, including Web animations that show the dynamics of the singularities in the complex plane, are presented.
Numerical Integration of Periodic Functions: A Few Examples
American Mathematical Monthly, Vol. 109, pp. 21--36 (2002)
It comes as a surprise to many students that the very basic trapezoidal and midpoint rules for numerical integration can yield more accurate results than more sophisticated methods such as the Simpson and Gauss-Legendre rules. This happens in particular when the integrand is smooth and periodic. In a series of examples explicit error formulas for trapezoidal/midpoint approximations to a few integrals are derived. These examples were chosen to demonstrate a variety of convergence behaviors, most of them quicker than what the standard error estimates predict.
Quadrature Rules Based on Partial Fraction Expansions (with D.P. Laurie)
Numerical Algorithms, Vol. 24, pp. 159--178 (2000) (online)
Quadrature rules are typically derived by requiring that all polynomials of a certain degree be integrated exactly. The nonstandard issue discussed here is the requirement that, in addition to the polynomials, the rule also integrates a set of prescribed rational functions exactly. Recurrence formulas for computing such quadrature rules are derived. In addition, Fejer's first rule, which is based on polynomial interpolation at Chebyshev nodes, is extended to integrate also rational functions with pre-assigned poles exactly. Numerical results, showing a favorable comparison with similar rules that have been proposed in the literature, are presented. An error analysis of a representative test problem is given.
A MATLAB Differentiation Matrix Suite (with S.C. Reddy)
ACM Transactions of Mathematical Software, Vol. 26, pp. 465--519 (2000)
A software suite consisting of seventeen Matlab functions for solving differential equations by the spectral collocation (a.k.a. pseudospectral) method is presented. It includes functions for computing derivatives of arbitrary order corresponding to Chebyshev, Hermite, Laguerre, Fourier, and sinc interpolants. Auxiliary functions are included for incorporating boundary conditions, performing interpolation using barycentric formulas, and computing roots of orthogonal polynomials. It is demonstrated how to use the package for solving eigenvalue, boundary value, and initial value problems arising in the fields of special functions, quantum mechanics, nonlinear waves, and hydrodynamic stability.
Spectral Methods Based on Non-Classical Orthogonal Polynomials
International Series on Numerical Mathematics, Vol. 131, pp. 238--251 (1999)
Spectral methods for solving differential equations of boundary value type have traditionally been based on classical orthogonal polynomials such as the Chebyshev, Legendre, Laguerre, and Hermite polynomials. In this numerical study we show that methods based on non-classical orthogonal polynomials may sometimes be more accurate. Examples include the solution of a two-point boundary value problem with a steep boundary layer and two Sturm-Liouville problems.
Algorithms for Parameter Selection in the Weeks Method for Inverting the Laplace Transform
SIAM Journal on Scientific Computing, Vol. 21, pp. 111-128 (1999)
The Weeks method is one of the most efficient numerical techniques for inverting the Laplace transform, provided two free parameters in the Laguerre expansion on which the method is based are chosen judiciously. However, there appears to be no theoretical estimates nor numerical algorithms for computing these parameters. The two algorithms presented in this paper fill that gap. Both algorithms aim to minimize a theoretical error estimate, and both are implemented with the FFT. The first algorithm is restricted to a certain class of transforms and the user is expected to provide information on the singularities of the transform. The second algorithm, though more expensive, is completely general---no singularity information is required. In challenging numerical tests both algorithms successfully predicted values of the parameters that are near optimal.
Comment on: ``A note on an integrable discretization of the nonlinear Schrödinger equation'' (with W. Black, B.M. Herbst)
Inverse Problems, Vol. 15, pp. 807-810 (1999)
In a recent paper [Inverse Problems, Vol. 13, (1997), 1121--1136], the author suggested a new implementation of an integrable discretization of the nonlinear Schrödinger equation due to Ablowitz and Ladik. The new implementation overcomes the non-locality of the original scheme. In this note it is shown, however, that the proposed scheme suffers from a severe grid-size restriction, which makes it unattractive from a practical viewpoint.
Finite Difference Methods for an AKNS Eigenproblem (with B.M. Herbst).
Mathematics and Computers in Simulation, Vol. 43, pp. 77-88 (1997)
We consider the numerical solution of the AKNS eigenproblem associated with the nonlinear Schrödinger equation. Four finite difference methods are considered: two standard schemes (forward and central differences), a discretization introduced by Ablowitz and Ladik, and a modified version of the latter scheme. By comparing these methods both numerically and theoretically we show that the modified Ablowitz-Ladik scheme has several desirable features. This includes the property that with a given number of gridpoints it approximates much larger sections of the spectrum than its rivals.
Computing the Hilbert Transform on the Real Line
Mathematics of Computation, Vol. 64, pp. 745-762 (1995)
We introduce a new method for computing the Hilbert transform on the real line. It is a collocation method, based on an expansion in rational eigenfunctions of the Hilbert transform operator, and implemented through the Fast Fourier Transform. An error analysis is given, and convergence rates for some simple classes of functions are established. Numerical tests indicate that the method compares favorably with existing methods.
Computation of the Complex Error Function
SIAM Journal on Numerical Analysis, Vol. 31, pp. 1497-1518 (1994) (A production error involving some figures was corrected in SIAM Journal on Numerical Analysis, Vol. 32, pp. 330-331 (1995))
We present rational expansions for computing the complex error function w(z) = exp(-z^2) erfc(-iz). These expansions have the following attractive properties: (1) They can be evaluated using a polynomial evaluation routine such as Horner's method, (2) the polynomial coefficients can be computed once and for all by a single FFT, and (3) high accuracy is achieved uniformly in the complex plane with only a small number of terms. Comparisons reveal that in some parts of the complex plane certain competitors may be more efficient. However, the difference in efficiency is never great, and our new algorithms are simpler than existing ones: a complete program takes eight lines of Matlab code.
Computing Integrals of the Complex Error Function
Proceedings of Symposia in Applied Mathematics, Vol. 48, pp. 403-407 (1994)
We present an efficient expansion for computing the integral of the complex error function. Special cases include the integrals of the Dawson and complementary error functions.
Theory and Applications of an Orthogonal Rational Basis Set
Proceedings of the Twentieth South African Symposium on Numerical Mathematics, Unhlanga Rocks, 4-6 July 1994.
This paper reviews the theory and applications of the functions
phi_n(x) = (1+ix)^n/(1-ix)^(n+1), which form a complete and orthogonal
basis set for L_2(R). It is shown how these functions may be used
for the computation of integral transforms and certain special functions.
The numerical solution of some differential equations is also discussed.
(scanned copy)
Numerical Simulation of Solitons and Dromions in the Davey-Stewartson System (with P.W. White)
Mathematics and Computers in Simulation, Vol. 37, pp. 469-479 (1994)
We extend the well-known split-step Fourier method for solving the nonlinear Schr\"odinger equation to the Davey-Stewartson system. We discuss both the elliptic-hyperbolic and hyperbolic-elliptic cases, as well as the implementation of boundary conditions. Numerical tests, including the simulation of rational solitons and dromions, are presented.
The eigenvalues of Hermite and rational spectral differentiation matrices
Numerische Mathematik, Vol. 61, pp. 409--431 (1992)
We derive expressions for the eigenvalues of spectral differentiation matrices for unbounded domains. In particular, we consider Galerkin and collocation methods based on Hermite functions as well as rational functions (a Fourier series combined with a cotangent mapping). We show that (i) first derivative matrices have purely imaginary eigenvalues and second derivative matrices have real and negative eigenvalues, (ii) for the Hermite method the eigenvalues are determined by the roots of the Hermite polynomials and for the rational method they are determined by the Laguerre polynomials, and (iii) the Hermite method has attractive stability properties in the sense of small condition numbers and spectral radii.
Pseudospectral Methods for the Benjamin-Ono Equation (with R.L. James)
Advances in Computer Methods for Partial Differential Equations VII, pp. 371-377 (1992)
We propose two pseudospectral methods for solving the Benjamin-Ono
equation. One is based on Fourier series, the other on rational
functions. We discuss various issues such as the implementation
of these methods and their conservation properties, and we
compare them numerically.
(scanned copy)
An adaptive algorithm for spectral computations on unbounded intervals (with A. Cloot)
Journal of Compuational Physics, Vol. 102, pp. 398-406 (1992)
We consider spectral methods for solving differential equations on unbounded domains. In particular, we consider a spectral rational function method and a method based on domain truncation combined with a Fourier series. Both methods contain a free parameter which determines the accuracy of the spectral method. When solving evolution equations the optimal parameter may change with time resulting in a significant loss of accuracy if the parameter is kept fixed. We develop an automatic algorithm which adapts the parameter at regular time levels to maintain optimum accuracy. In numerical tests we sound that the proposed method is robust and efficient.
Dynamics of Semi-Discretizations of the Defocusing Nonlinear Schrödinger Equation (with B.M. Herbst and M.J. Ablowitz)
IMA Journal of Numerical Analysis, Vol. 11, pp.~539--552 (1991)
It has been demonstrated that the nonlinear Schrödinger equation is sensitive to discretizations. In the focussing case this is due to the homoclinic structure associated with the NLS equation. In this paper we show that various numerical schemes for the defocusing case are also prone to instabilities, although not as severe as those of the focusing equation. An integrable discretization due to Ablowitz and Ladik does not suffer from the same instabilities. However, it is shown that it develops a focusing singularity if a threshold condition is exceeded. Numerical examples illustrating the phenomena pertaining to the defocsing equation are given.
Two results on polynomial interpolation in equally spaced points (with L.N. Trefethen)
Journal of Approximation Theory, Vol. 65, pp. 247-260 (1991)
We present two results that quantify the poor behavior of polynomial interpolation in n equally spaced points. First, in band-limited interpolation of complex exponential functions exp(i a x) (with a real), the error decreases to 0 as n--> infinity if and only if a is small enough to provide at least 6 points per wavelength. Second, the Lebesgue constant L_n (supremum norm of the n-th degree interpolation operator) satisfies lim_{n --> \infty} L_n^(1/n) = 2. Both of these results are more than fifty years old, but they are generally unknown to approximation theorists.
Spectral methods and mappings for evolution equations on the infinite line (with A. Cloot)
Computer Methods in Applied Mechanics and Engineering, Vol. 80, pp. 467-481 (1990).
Also published in the Proceedings of ICOSAHOM '89, Como, Italy, 26-29 June 1989.
Mapping methods for improving the accuracy of Fourier psuedospectral differentiation on the infinite line are considered. For a certain class of functions, optimal mapping parameters are obtained by balancing the error due to boundary effects with the internal resolution error. Estimates for the maximal rate of convergence yielded by each mapping are also derived. In the case of evolution equations, the effect of these mappings on the stability of the time integration is investigated. The robustness of the optimal mapping parameters with respect to changes in the solution during the time evolution is examined numerically.
A numerical study of the nonlinear Schrödinger equation involving quintic terms (with A. Cloot and B.M. Herbst)
Journal of Computational Physics, Vol. 86, pp. 127-146 (1990)
The cubi-quintic Schrödinger equation is known to possess solutions that grow unboundedly in finite time. By exploring its conservation properties we derive sufficient conditions for bounded solutions. The computation of solutions near the critical threshold poses difficulties, since the number of active Fourier-components increase dramatically, resulting in steep temporal and spatial gradients. To overcome this difficulty we propose an efficient pseudospectral scheme which adaptively adjust the number of degrees of freedom.
The eigenvalues of second-order spectral differentiation matrices (with L.N. Trefethen)
SIAM Journal on Numerical Analysis, Vol. 25, pp. 1279-1298 (1988)
The eigenvalues of the pseudospectral second derivative matrix with homogeneous Dirchlet boundary conditions are important in many applications of spectral methods. This paper investigates some of their properties. Numerical results show that a certain fraction of the eigenvalues approximate the continuous operator very accurately, but the errors in the remaining ones are large. It is demonstrated that the inaccurate eigenvalues correspond to those eigenfunctions of the continuous operator that cannot be resolved by polynomial interpolation in the spectral grid. In particular, it is proved that pi points on average per wavelength are sufficient for successful interpolation of the eigenfunctions of the continuous operator in a Chebyshev distribution of nodes, and six points per wavelength for a uniform distribution. These results are in agreement with the observed fractions of accurate eigenvalues. By using the characteristic polynomial, a bound on the spectral radius of the differentiation matrix is derived that is accurate to 2% or better. The effect of filtering on the eigenvalues is examined numerically.
Recurrence in semi-discrete approximations of the nonlinear Schrödinger equation (with B.M. Herbst)
SIAM Journal of Scientific and Statistical Computing, Vol. 8, pp. 988-1004 (1987)
Recurrence in the nonlinear Schrödinger equation is studied by means of three semidiscrete approximations, namely a spectral, a pseudospectral, and a finite difference method. These methods are analyzed for conservation and stability properties. A linear stability analysis is presented, as well as a nonlinear two-mode analysis. A comparison of the results reveals some disadvantages of the pseudospectral and finite difference methods. The role of aliasing is examined and some numerical results are given.
Split-step methods for the solution of the nonlinear Schrödinger equation (with B.M. Herbst)
SIAM Journal on Numerical Analysis, Vol. 23, pp. 485-507 (1986)
A split-step method is used to discretize the time variable for the numerical solution of the nonlinear Schrödinger equation. The space variable is discretized by means of a finite difference and a Fourier method. These methods are analyzed with respect to various physical properties represented in the NLS. In particular, it is shown how a conservation law, dispersion, and stability are reflected in the numerical schemes. Analytical and numerical instabilities of wave train solutions are identified and conditions which avoid the latter are derived. These results are demonstrated by numerical examples.
On the stability of the nonlinear Schrödinger equation (with B.M. Herbst, A.R. Mitchell)
Journal of Computational Physics, Vol. 60, pp. 263-281 (1985)
The analytical and numerical stability properties of the nonlinear Schrödinger equation are investigated. It is shown how the analytical instabilities are reflected in the semi- and fully discrete numerical schemes. Numerical results of nonlinear blowup and recurrence are given.
Last updated: January 10, 2020