Parabolic and Hyperbolic Contours for Computing the Bromwich Integral

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.\ud \ud JACW was supported by the National Research Foundation in South Africa under grant FA200503230001

Computing the Gamma function using contour integrals and rational approximations

Some of the best methods for computing the gamma function are based on numerical evaluation of Hankel's contour integral. For example, Temme evaluates this integral based on steepest-decent contours by the trapezoid rule. Here we investigate a different approach to the integral: the application of the trapezoid rule on Talbot-type contours using optimal parameters recently derived by Weideman for computing inverse Laplace transforms. Relatedly, we also investigate quadrature formulas derived from best approximations to exp(z) on the negative real axis, following Cody, Meinardus and Varga. The two methods are closely related and both converge geometrically. We find that the new methods are competitive with existing ones, even though they are based on generic tools rather than on specific analysis of the gamma function

Numerical evaluation of two and three parameter Mittag-Leffler functions

The Mittag-Leffler (ML) function plays a fundamental role in fractional calculus but very few methods are available for its numerical evaluation. In this work we present a method for the efficient computation of the ML function based on the numerical inversion of its Laplace transform (LT): an optimal parabolic contour is selected on the basis of the distance and the strength of the singularities of the LT, with the aim of minimizing the computational effort and reduce the propagation of errors. Numerical experiments are presented to show accuracy and efficiency of the proposed approach. The application to the three parameter ML (also known as Prabhakar) function is also presented.Comment: Accepted for publication in SIAM Journal on Numerical Analysi

A robust spectral method for solving Heston’s model

In this paper, we consider the Heston’s volatility model (Heston in Rev. Financ. Stud. 6: 327–343, 1993]. We simulate this model using a combination of the spectral collocation method and the Laplace transforms method. To approximate the two dimensional PDE, we construct a grid which is the tensor product of the two grids, each of which is based on the Chebyshev points in the two spacial directions. The resulting semi-discrete problem is then solved by applying the Laplace transform method based on Talbot’s idea of deformation of the contour integral (Talbot in IMA J. Appl. Math. 23(1): 97–120, 1979)

Contour integral method for European options with jumps

We develop an efficient method for pricing European options with jump on a single asset. Our approach is based on the combination of two powerful numerical methods, the spectral domain decomposition method and the Laplace transform method. The domain decomposition method divides the original domain into sub-domains where the solution is approximated by using piecewise high order rational interpolants on a Chebyshev grid points. This set of points are suitable for the approximation of the convolution integral using Gauss–Legendre quadrature method. The resulting discrete problem is solved by the numerical inverse Laplace transform using the Bromwich contour integral approach. Through rigorous error analysis, we determine the optimal contour on which the integral is evaluated. The numerical results obtained are compared with those obtained from conventional methods such as Crank–Nicholson and finite difference. The new approach exhibits spectrally accurate results for the evaluation of options and associated Greeks. The proposed method is very efficient in the sense that we can achieve higher order accuracy on a coarse grid, whereas traditional methods would required significantly more time-steps and large number of grid points.Web of Scienc

Evaluating matrix functions for exponential integrators via Carathéodory-Fejér approximation and contour integrals

Among the fastest methods for solving stiff PDE are exponential integrators, which require the evaluation of $f(A)$, where $A$ is a negative definite matrix and $f$ is the exponential function or one of the related ``$\varphi$ functions'' such as $\varphi_1(z) = (e^z-1)/z$. Building on previous work by Trefethen and Gutknecht, Gonchar and Rakhmanov, and Lu, we propose two methods for the fast evaluation of $f(A)$ that are especially useful when shifted systems $(A+zI)x=b$ can be solved efficiently, e.g. by a sparse direct solver. The first method method is based on best rational approximations to $f$ on the negative real axis computed via the Carathéodory-Fejér procedure, and we conjecture that the accuracy scales as $(9.28903\dots)^{-2n}$, where $n$ is the number of complex matrix solves. In particular, three matrix solves suffice to evaluate $f(A)$ to approximately six digits of accuracy. The second method is an application of the trapezoid rule on a Talbot-type contour

The exponentially convergent trapezoidal rule

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