Today, symbolic mathematical computation is taken for granted as part of the scientific infrastructure, but it has not always been so. This paper provides a historical survey of the discipline’s formative decade, 1965–1975, viewed from a 50 year perspective. This span of years saw the evolution from a few specialized programs with naive algorithms to integrated systems with substantial capabilities. We highlight some of the important early figures in the field and the innovations upon which the current generation of systems and algorithms are built. By revisiting a period unfamiliar to most current readers, this survey aims to shed light on once-pressing issues that are now largely resolved and to highlight how some of today’s challenges were recognized earlier than expected.

Given the set of discrete solution points or nodes, called the skeleton, generated by an ODE solver, we study the problem of fitting a curve passing through the nodes in the skeleton minimizing a norm of the residual vector of the ODE. We reformulate this interpolation problem as a multi-stage optimal control problem and, for the minimization of two different norms, we apply the associated maximum principle to obtain the necessary conditions of optimality. We solve the problem analytically for the Dahlquist test problem and a variant of the leaky bucket problem, in terms of the given skeleton. We also consider the Van der Pol equation, for which we obtain interpolating curves with minimal residual norms by numerically solving a direct discretization of the problem through optimization software. With the skeletons obtained by various ODE solvers of MATLAB, we make comparisons between the residuals obtained by our approach and those obtained by the MATLAB function deval.

The classical WKB method for solving singularly perturbed linear differential equations has never, as far as we know, been looked at from the structured backward error analysis point of view. This is somewhat surprising because a simple computation shows that for some important problems, the WKB method gives the exact solution of a problem of the same structure that can be expressed in finitely many terms. This kind of analysis can be extremely useful in assessing the validity of a solution provided by the WKB method. In this chapter, we show how to do this and explore some of the consequences, which include a new iterative algorithm to improve the quality of the WKB solution. We also explore a new hybrid method where the potential is approximated by Chebyshev polynomials, which can be implemented in a few lines of Chebfun.

The definition of Stirling numbers for complex arguments proposed by Flajolet Prodinger is studied. Stirling numbers are defined by contour integrals in the complex plane. In addition to symbolic special cases, methods for evaluating the integrals numerically are presented. Numerical integration methods based on exponentially convergent trapezoidal schemes are implemented; optimizations of the calculations are explored.

The book under review was originally published under the auspices of the National Research Council in 1933 (the year John was born), and it was republished as a Dover edition in 1956 (three years before Rob was born). At 108 pages—including title page, preface, table of contents, and index—it’s very short. Even so, it contains a significant amount of information that was of technical importance for its time and is of historical importance now.

Stirling cycle numbers and Stirling partition numbers have many combinatorial applications. A symmetric matrix built from Stirling cycle numbers, which is known to be totally nonnegative, appears as https://oeis.org/A137854 in the Online Encyclopedia of Integer Sequences. In this paper we give analytical bidiagonal factorings of these matrices.

The Hermite interpolation problem—defined in the article text—is more complicated than the Lagrange interpolation problem—also defined there—and occurs less frequently in practice. But it does occur, and solving it is occasionally useful. Solutions have been reinvented many times since the problem was first posed and solved in 1878 by Charles Hermite. This article shows how the barycentric forms of the solution, invented about a hundred years after Hermite, work. All one needs to do is to compute a partial fraction expansion by a numerically stable method, and this gives us numerically stable and efficient forms to evaluate the Hermite interpolational polynomial. I describe the Maple program BHIP and some of its ancillary routines (available for download from the Maple Cloud, by clicking on the link below the link to the article PDF, to the right of this abstract), and mention the equivalent Matlab versions genbarywts and hermiteeval. I also compare to some less numerically stable and less efficient approaches. I also show how to find the roots of a Hermite interpolational polynomial by constructing a companion matrix pencil with the routine CMP, which does not change the polynomial basis, and then using standard software to compute the generalized eigenvalues of the pencil, which then give us the roots of the polynomial.

The numerical solution of initial-value problems (IVP) for ordinary differential equations (ODE) is at this time a mature subject, with many high-quality codes freely available. Second-order linear equations without singularities are an especially simple class of problems to solve, even more so if only a single scalar equation such as the Mathieu equation y” + (a - 2qcos 2x)y = 0 is being considered. Nonetheless, the topic is not yet exhausted, and this paper considers the case of writing an efficient arbitrary-precision code for the solution of such equations. For this purpose, an implicit Hermite–Obreshkov method attains nearly spectral accuracy at a cost only polynomial in the number of bits of accuracy requested. This is interesting for the Mathieu equation in particular because the solutions can be highly oscillatory of variable frequency and be highly ill-conditioned. This paper reports on the details of the prototype Maple implementation of the method and summarizes the approximation theoretic results justifying the choice of a balanced Hermite–Obreshkov method including its backward stability and decent Lebesgue constants. This method may be of especial interest for the solution of so-called D-finite equations, for which Taylor series coefficients up to degree m are available at cost only O ( m ), instead of the more usual O(m^2) . This paper celebrates the happy occasion of the 90th birthday of John C. Butcher.

A “blendstring” is a piecewise polynomial interpolant with high-degree two-point Hermite interpolational polynomials on each piece, analogous to a cubic spline. Blendstrings are smoother and can be more accurate than cubic splines, and can be used to represent smooth functions on a line segment or polygonal path in the complex plane. I sketch some properties of blendstrings, including efficient methods for evaluation, differentiation, and integration, as well as a prototype Maple implementation. Blendstrings can be differentiated and integrated exactly and can be combined algebraically. I also show applications of blendstrings to solving differential equations and computing Mathieu functions and generalized Mathieu eigenfunctions.

This paper discusses some of the philosophical and historical underpinnings of the talk “The Mathieu Functions: Computational and Historical Perspectives” given at the Maple Conference 2022. I also comment on the role Problem Solving Environments (PSEs) play in curating the computational knowledge of the 19th century, which is so necessary for us to think of special functions as answers rather than questions.

Abstract The numerical solution of initial-value problems (IVP) for ordinary differentialequations (ODE) is at this time a mature subject, with many high-quality codes freely available.Second-order linear equations without singularities are an especially simple class ofproblems to solve, even more so if only a single scalar equation such as the Mathieu equation y′′+(a−2qcos2x)y=0 is being considered. Nonetheless, the topic is not yet exhausted, andthis paper considers the case of writing an efficient arbitrary-precision code for the solutionof such equations. For this purpose, an implicit Hermite–Obreshkov method attains nearlyspectral accuracy at at cost only polynomial in the number of bits of accuracy requested.This is interesting for the Mathieu equation in particular because the solutions can be highlyoscillatory of variable frequency and be highly ill-conditioned. This paper reports on the details of the prototype Maple implementation of the method,and summarizes the approximation theoretic results justifying the choice of a balancedHermite–Obreshkov method including its backward stability and decent Lebesgue constants.This method may be of especial interest for the solution of so-called D-finite equations,for which Taylor series coefficients up to degree m are available at cost only O(m),instead of the more usual O(m2). This paper celebrates the happy occasion of the 90th birthday of John C. Butcher. Mathematics Subject Classification (2000) 65L04 · 33F05 · 65D15

The authors have been using a largely algebraic form of ``computational discovery'' in various undergraduate classes at their respective institutions for some decades now to teach pure mathematics, applied mathematics, and computational mathematics. This paper describes what we mean by ``computational discovery,'' what good it does for the students, and some specific techniques that we used.

The recursively-constructed family of Mandelbrot matrices M-n for n = 1, 2,... have nonnegative entries (indeed just 0 and 1, so each M-n can be called a binary matrix) and have eigenvalues whose negatives -lambda = c give periodic orbits under the Mandelbrot iteration, namely z(k) = z(k-1)(2) + c with z(0) = 0, and are thus contained in the Mandelbrot set. By the Perron-Frobenius theorem, the matrices M-n have a dominant real positive eigenvalue, which we call rho(n). This article examines the eigenvector belonging to that dominant eigenvalue and its fractal-like structure, and similarly examines (with less success) the dominant singular vectors of M-n from the singular value decomposition.

In this paper we introduce an apparently new spline-like interpolant that I call a compact cubic interpolant or compact cubic spline; this is similar to a cubic spline introduced in 1972 by Swartz and Varga, but has higher order accuracy at the edges. We argue that for nearly uniform meshes the compact cubic approach offers some potential advantages, and offers a simple way to treat the edge conditions, relieving the user of the burden of deciding to use one of the three standard options: free (natural), complete (clamped), or “not-a-knot” conditions. Finally, we establish that the matrices defining the compact cubic splines (equivalently, the fourth-order compact finite difference formulæ) are totally nonnegative, if all mesh widths are the same sign, for instance if the mesh is real and nodes are numbered in increasing order.

The authors have been using a largely algebraic form of “computational discovery” in various undergraduate classes at their respective institutions for some decades now to teach pure mathematics, applied mathematics, and computational mathematics. This paper describes what we mean by “computational discovery,” what good it does for the students, and some specific techniques that we used.

One useful standard method to compute eigenvalues of matrix polynomials P(z)∈ C n x n [z] of degree at most ℓ in z (denoted of grade ℓ, for short) is to first transform P(z) to an equivalent linear matrix polynomial L(z)=zB-A, called a companion pencil, where A and B are usually of larger dimension than P(z) but L(z) is now only of grade 1 in z. The eigenvalues and eigenvectors of L(z) can be computed numerically by, for instance, the QZ algorithm. The eigenvectors of P(z), including those for infinite eigenvalues, can also be recovered from eigenvectors of L(z) if L(z) is what is called a "strong linearization'' of P(z). In this paper we show how to use algorithms for computing the Hermite Normal Form of a companion matrix for a scalar polynomial to direct the discovery of unimodular matrix polynomial cofactors E(z) and F(z) which, via the equation E(z)L(z)F(z) = diag(P(z), In, …, I_n), explicitly show the equivalence of P(z) and P(z). By this method we give new explicit constructions for several linearizations using different polynomial bases. We contrast these new unimodular pairs with those constructed by strict equivalence, some of which are also new to this paper. We discuss the limitations of this experimental, computational discovery method of finding unimodular cofactors.

Differentiation matrices are in wide use in numerical algorithms, although usually studied in an ad hoc manner. We collect here in this review paper elementary properties of differentiation matrices for univariate polynomials expressed in various bases, including orthogonal polynomial bases and non-degree-graded bases such as Bernstein bases and Lagrange and Hermite interpolational bases. We give new explicit formulations, and new explicit formulations for the pseudo-inverses which help to understand antidifferentiation, of many of these matrices. We also give the unique Jordan form for these (nilpotent) matrices and a new unified formula for the transformation matrix.

We consider the problem of symbolic-numeric integration of symbolic functions, focusing on rational functions. Using a hybrid method allows the reliable yet efficient computation of symbolic antiderivatives while avoiding issues of ill-conditioning to which numerical methods are susceptible. We propose two alternative methods for exact input that compute the rational part of the integral using Hermite reduction and then compute the transcendental part two different ways using a combination of exact integration and efficient numerical computation of roots. The symbolic computation is done within bpas, or Basic Polynomial Algebra Subprograms, which is a highly optimized environment for polynomial computation on parallel architectures, while the numerical computation is done using the highly optimized multiprecision rootfinding package MPSolve. We provide for both algorithms computable expressions for the first-order term of a structured forward and backward error and show how, away from singularities, tolerance proportionality is achieved by adjusting the precision of the rootfinding tasks.

We look at two classical examples in the theory of numerical analysis, namely the Runge example for interpolation and Wilkinson's example (actually two examples) for rootfinding. We use the modern theory of backward error analysis and conditioning, as instigated and popularized by Wilkinson, but refined by Farouki and Rajan. By this means, we arrive at a satisfactory explanation of the puzzling phenomena encountered by students when they try to fit polynomials to numerical data, or when they try to use numerical rootfinding to find polynomial zeros. Computer algebra, with its controlled, arbitrary precision, plays an important didactic role.