A New Algorithm for Computing the Actions of Trigonometric and Hyperbolic Matrix Functions

A new algorithm is derived for computing the actions $f(tA)B$ and $f(tA^{1/2})B$, where $f$ is cosine, sinc, sine, hyperbolic cosine, hyperbolic sinc, or hyperbolic sine function. $A$ is an $n\times n$ matrix and $B$ is $n\times n_0$ with $n_0 \ll n$. $A^{1/2}$ denotes any matrix square root of $A$ and it is never required to be computed. The algorithm offers six independent output options given $t$, $A$, $B$, and a tolerance. For each option, actions of a pair of trigonometric or hyperbolic matrix functions are simultaneously computed. The algorithm scales the matrix $A$ down by a positive integer $s$, approximates $f(s^{-1}tA)B$ by a truncated Taylor series, and finally uses the recurrences of the Chebyshev polynomials of the first and second kind to recover $f(tA)B$. The selection of the scaling parameter and the degree of Taylor polynomial are based on a forward error analysis and a sequence of the form $\|A^k\|^{1/k}$ in such a way the overall computational cost of the algorithm is optimized. Shifting is used where applicable as a preprocessing step to reduce the scaling parameter. The algorithm works for any matrix $A$ and its computational cost is dominated by the formation of products of $A$ with $n\times n_0$ matrices that could take advantage of the implementation of level-3 BLAS. Our numerical experiments show that the new algorithm behaves in a forward stable fashion and in most problems outperforms the existing algorithms in terms of CPU time, computational cost, and accuracy.Comment: 4 figures, 16 page

A DPG-based time-marching scheme for linear hyperbolic problems

The Discontinuous Petrov-Galerkin (DPG) method is a widely employed discretization method for Partial Di fferential Equations (PDEs). In a recent work, we applied the DPG method with optimal test functions for the time integration of transient parabolic PDEs. We showed that the resulting DPG-based time-marching scheme is equivalent to exponential integrators for the trace variables. In this work, we extend the aforementioned method to time-dependent hyperbolic PDEs. For that, we reduce the second order system in time to first order and we calculate the optimal testing analytically. We also relate our method with exponential integrators of Gautschi-type. Finally, we validate our method for 1D/2D + time linear wave equation after semidiscretization in space with a standard Bubnov-Galerkin method. The presented DPG-based time integrator provides expressions for the solution in the element interiors in addition to those on the traces. This allows to design di fferent error estimators to perform adaptivity

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

Accurate Approximation of the Matrix Hyperbolic Cosine Using Bernoulli Polynomials

[EN] This paper presents three different alternatives to evaluate the matrix hyperbolic cosine using Bernoulli matrix polynomials, comparing them from the point of view of accuracy and computational complexity. The first two alternatives are derived from two different Bernoulli series expansions of the matrix hyperbolic cosine, while the third one is based on the approximation of the matrix exponential by means of Bernoulli matrix polynomials. We carry out an analysis of the absolute and relative forward errors incurred in the approximations, deriving corresponding suitable values for the matrix polynomial degree and the scaling factor to be used. Finally, we use a comprehensive matrix testbed to perform a thorough comparison of the alternative approximations, also taking into account other current state-of-the-art approaches. The most accurate and efficient options are identified as results.This research was supported by the Vicerrectorado de Investigacion de la Universitat Politecnica de Valencia (PAID-11-21).Alonso Abalos, JM.; Ibáñez González, JJ.; Defez Candel, E.; Alvarruiz Bermejo, F. (2023). Accurate Approximation of the Matrix Hyperbolic Cosine Using Bernoulli Polynomials. Mathematics. 11(3):1-22. https://doi.org/10.3390/math1103052012211

On computing high-dimensional Riemann theta functions

Riemann theta functions play a crucial role in the field of nonlinear Fourier analysis, where they are used to realize inverse nonlinear Fourier transforms for periodic signals. The practical applicability of this approach has however been limited since Riemann theta functions are multi-dimensional Fourier series whose computation suffers from the curse of dimensionality. In this paper, we investigate several new approaches to compute Riemann theta functions with the goal of unlocking their practical potential. Our first contributions are novel theoretical lower and upper bounds on the series truncation error. These bounds allow us to rule out several of the existing approaches for the high-dimension regime. We then propose to consider low-rank tensor and hyperbolic cross based techniques. We first examine a tensor-train based algorithm which utilizes the popular scaling and squaring approach. We show theoretically that this approach cannot break the curse of dimensionality. Finally, we investigate two other tensor-train based methods numerically and compare them to hyperbolic cross based methods. Using finite-genus solutions of the Korteweg–de Vries (KdV) and nonlinear Schrödinger equation (NLS) equations, we demonstrate the accuracy of the proposed algorithms. The tensor-train based algorithms are shown to work well for low genus solutions with real arguments but are limited by memory for higher genera. The hyperbolic cross based algorithm also achieves high accuracy for low genus solutions. Its novelty is the ability to feasibly compute moderately accurate solutions (a relative error of magnitude 0.01) for high dimensions (up to 60). It therefore enables the computation of complex inverse nonlinear Fourier transforms that were so far out of reach

Spectral methods for CFD

One of the objectives of these notes is to provide a basic introduction to spectral methods with a particular emphasis on applications to computational fluid dynamics. Another objective is to summarize some of the most important developments in spectral methods in the last two years. The fundamentals of spectral methods for simple problems will be covered in depth, and the essential elements of several fluid dynamical applications will be sketched

Computing the Wave-Kernel Matrix Functions

We derive an algorithm for computing the wave-kernel functions cosh \surd A and sinhc\surd A for an arbitrary square matrix A, where sinhcz = sinh(z)/z. The algorithm is based on Pad\'e approximation and the use of double angle formulas. We show that the backward error of any approximation to cosh \surd A can be explicitly expressed in terms of a hypergeometric function. To bound the backward error we derive and exploit a new bound for \| Ak\| 1/k that is sharper than one previously obtained by Al-Mohy and Higham [SIAM J. Matrix Anal. Appl., 31 (2009), pp. 970-- 989]. The amount of scaling and the degree of the Pad\'e approximant are chosen to minimize the computational cost subject to achieving backward stability for cosh \surd A in exact arithmetic. Numerical experiments show that the algorithm behaves in a forward stable manner in floating-point arithmetic and is superior in this respect to the general purpose Schur--Parlett algorithm applied to these functions

Secular dynamics of a planar model of the Sun-Jupiter-Saturn-Uranus system; effective stability into the light of Kolmogorov and Nekhoroshev theories

We investigate the long-time stability of the Sun-Jupiter-Saturn-Uranus system by considering a planar secular model, that can be regarded as a major refinement of the approach first introduced by Lagrange. Indeed, concerning the planetary orbital revolutions, we improve the classical circular approximation by replacing it with a solution that is invariant up to order two in the masses; therefore, we investigate the stability of the secular system for rather small values of the eccentricities. First, we explicitly construct a Kolmogorov normal form, so as to find an invariant KAM torus which approximates very well the secular orbits. Finally, we adapt the approach that is at basis of the analytic part of the Nekhoroshev's theorem, so as to show that there is a neighborhood of that torus for which the estimated stability time is larger than the lifetime of the Solar System. The size of such a neighborhood, compared with the uncertainties of the astronomical observations, is about ten times smaller.Comment: 31 pages, 2 figures. arXiv admin note: text overlap with arXiv:1010.260