Quadrature Rules Based on Partial Fraction Expansions (with D.P. Laurie)
To appear in Numerical Algorithms
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)
To appear in ACM Transactions of Mathematical Software
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 Schrodinger 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 Schrodinger 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\"odinger 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.
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.
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\"odinger 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\"odinger 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\"odinger 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\"odinger 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\"odinger equation (with B.M. Herbst)
SIAM Journal of Scientific and Statistical Computing, Vol. 8, pp. 988-1004 (1987)
Recurrence in the nonlinear Schr\"odinger 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\"odinger 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\"odinger 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\"odinger 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\"odinger 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: July 21, 1999