The paper exhibits several standing modes of a 5:3 elliptic drum computed using Mathieu functions. To match the boundary conditions, I used Newton's method on the appropriate modified Mathieu equation using the Squire-Trapp formula for computing derivatives. I tabulate the requisite values of the parameter q for these low-frequency modes.

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.

Let A be a square matrix with real entries. The spread of A is defined as the maximum of the distances among the eigenvalues of A. Let S_m[a,b] denote the set of all m× m symmetric matrices with entries in the real interval [a,b] and let S_m{a,b} be the subset of S_m[a,b] of Bohemian matrices with population from only the extremal elements {a,b}. S. M. Fallat and J. J. Xing in 2012 proposed the following conjecture: the maximum spread in S_m[a,b] is attained by a rank 2 matrix in S_m{a,b}. X. Zhan had proved previously that the conjecture was true for S_m[-a,a] with a>0. We will show how to interpret this problem geometrically, via polynomial resultants, in order to be able to treat this conjecture from a computational point of view. This will allow us to prove that this conjecture is true for several formerly open cases.

Previous work has made use of a parameterized plane curve polynomial representation for mathematical handwriting, with the polynomials represented in a Legendre or Legendre-Sobolev graded basis. This provides a compact geometric representation for the digital ink. Preliminary results have also been shown for Chebyshev and Chebyshev-Sobolev bases. This article explores the trade-offs between basis choice and polynomial degree to achieve accurate modeling with a low computational cost. To do this, we consider the condition number for polynomial evaluation in these bases and bound how the various inner products give norms for the variations between symbols.

This paper explores backward error analysis for numerical solutions of ordinary differential equations, particularly focusing on chaotic systems. Three approaches are examined: residual assessment, the method of modified equations, and shadowing. We investigate how these methods explain the success of numerical simulations in capturing the behavior of chaotic systems, even when facing issues like spurious chaos introduced by numerical methods or suppression of chaos by numerical methods. Finally, we point out an open problem, namely to explain why the statistics of long orbits are usually correct, even though we do not have a theoretical guarantee why this should be so.

Watson's lemma is a famous technique for approximate evaluation of integrals with a large parameter. It's easy to implement in Maple, but there's a surprise: Maple's implementation is actually stronger than the lemma, because int and series are so strong. In consequence, we went back and proved a stronger lemma than what had been known. The lemma can improve the speed of evaluating integrals, and we give an extended example with evalhf, the Maple compiler, and the Threads package.

The classical WKB method (also known as the WKBJ method, the LG method, or the phase integral method) for solving singularly perturbed linear differential equations has never, as far as we know, been looked at from the structured backward error (BEA) 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 paper 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.

An exponential tower is obtained by repeatedly raising a number to itself, as in xˣ, x^{x^x}, and so on. When formulated as an iterative process, it may converge to a single value, diverge to infinity, cycle between multiple values, or wander forever. Representing this property graphically leads to a beautiful fractal object. In order to plot the fractal, a method was devised to handle problems of numerical overflow and underflow that appear during the computation.

The infinite exponential tower is studied through the associated iteration c₁ = 0 and cₙ₊₁ = eᶜₙ λ, for complex λ. For a subset of λ values, the sequence displays stable 2-cycles, that is to say as n → ∞ we observe that the odd subsequence c₂ₙ₋₁ → A whereas the even subsequence c₂ₙ → B, with A ≠ B. Thus, A and B obey B=eᴬ λ and A = eᴮ λ. Numerical investigations of the 2-cycles use a further transformation ζexp(-ζ) = λ = ln(z) and the set of ζ values corresponding to 2-cycles has a curious shape, reminding us of pictures of insect larva; the region has sharply scalloped edges. This paper gives an analytic expression for the edges of the 2-cycle region and a complete explanation of the cusps on the boundary that give the scalloped look.

Some time in the early 2000's, I extended the routine CompanionMatrix in the LinearAlgebra package to compute what are called linearizations of what are called matrix polynomials. These are just univariate polynomials with matrix coefficients; isomorphically, these are matrices with univariate polynomial entries. Linearizations can be used to solve multivariate systems of equations by a number of techniques, which are `"well-known" in the sense that they are in books and papers. However well-known they are, they deserve to be better-known, and this expository paper gives some examples of some of the methods that can be used. Think of this as an extended help page for the code (which, if I am honest, is long overdue for an upgrade).

Arterial deformations arise in blood flow when surrounding tissue invades the space available for a blood vessel to maintain its circular cross section, the most immediate effects being a reduction in blood flow and redistribution of shear stress. Here we consider deformations from circular to elliptic cross sections. Solution of this problem in steady flow is fairly straightforward. The focus in the present paper is on pulsatile flow where the change from circular to elliptic cross sections is associated with a transition in the character of the equations governing the flow from Bessel to Mathieu equations. The main aim of our study is to examine the hemodynamic consequences of the change from circular to elliptic cross sections and on possible implications of this change in blood flow regulation. The study of this problem has been hampered in the past because of difficulties involved in the solution of the governing equations. In the present study we describe methods we have used to overcome some of these difficulties and present a comprehensive set of results based on these methods. In particular, vessel deformation is examined under two different conditions relevant to blood flow regulation: (i) keeping cross sectional area constant and (ii) keeping cross sectional circumference constant. The results provide an important context for the mechanism of neurovascular control of blood flow under the pathological conditions of vessel deformation. The difficulty which has characterized this problem is that it involves elements of mathematics which are well outside the scope of a clinical/physiological study, while it actually involves clinical/physiological elements which are well outside the scope of a mathematical study. We hope that the context which we provide in this paper helps resolve this difficulty.

This paper outlines our ideas on how to teach linear algebra in a mechanized mathematical environment, and discusses some of our reasons for thinking that this is a better way to teach linear algebra than the “old fashioned way”. We discuss some technological tools such as Maple, Matlab, Python, and Jupyter Notebooks, and some choices of topics that are especially suited to teaching with these tools. The discussion is informed by our experience over the past thirty or more years teaching at various levels, especially at the University of Western Ontario.

This note is a brief remembrance of our friend Eugenio Roanes–Lozano, who left us too soon.

The chaos game representation (CGR) is an interesting method to visualize one-dimensional sequences. In this paper, we show how to construct a chaos game representation. The applications mentioned here are biological, in which CGR was able to uncover patterns in DNA or proteins that were previously unknown. We also show how CGR might be introduced in the classroom, either in a modelling course or in a dynamical systems course. Some sequences that are tested are taken from the Online Encyclopedia of Integer Sequences, and others are taken from sequences that arose mainly from a course in experimental mathematics.

A "blend" is a two-point Hermite interpolational polynomial, typically of quite high degree. This note shows that implementing them in a double Horner evaluation scheme has good backward error, and also shows that the Lebesgue constant for a balanced blend or nearly balanced blend on the interval [0,1] is bounded by 2, independently of the grade or degree of the approximation. On [-1,1], which is a more natural interval for comparison, it is of course unbounded, but grows only like 2√(m/π) where 2m+1 is the grade of approximation. I also show that the quadrature schemes for balanced blends amplify errors only by O( ln(m) ).

The two most famous quotes about divergent series are Abel's "Divergent series are the invention of the devil, and it is shameful to base on them any demonstration whatsoever," and Heaviside's "This series is divergent, therefore we may be able to do something with it." Today a lot more is known about divergent series than in either's day, so we can say now that, on balance, Heaviside wins, and we now have plenty of license to use divergent series. This article talks about some "well-known" methods (that is, well-known to experts) to do so, and in particular talks about some of the devilishly good features of evalf/Sum, long one of my favourite tools in Maple. But Abel had a point, too, and we'll see some "shameful" things, which will give the reader some necessary caution to go along with their license.

We discuss the parallelization of algorithms for solving polynomial systems by way of triangular decomposition. The Triangularize algorithm proceeds through incremental intersections of polynomials to produce different components (points, curves, surfaces, etc.) of the solution set. Independent components imply the opportunity for concurrency. This “component-level” parallelization of triangular decompositions, our focus here, belongs to the class of dynamic irregular parallelism. Potential parallel speed-up depends only on geometrical properties of the solution set (number of components, their dimensions and degrees); these algorithms do not scale with the number of processors. To manage the irregularities of component-level parallelization we combine different concurrency patterns: map, workpile, producer-consumer, pipeline, and fork-join. We report on our implementation in the freely available BPAS library. Comprehensive experimentation with thousands of polynomial systems yields examples with up to 10.8× speed-up on a 12-core machine.

Over the past thirty years or so the authors have been teaching various programming for mathematics courses at our respective Universities, as well as incorporating computer algebra and numerical computation into traditional mathematics courses. These activities are, in some important ways, natural precursors to the use of Artificial Intelligence in Mathematics Education. This paper reflects on some of our course designs and experiences and is therefore a mix of theory and practice. Underlying both is a clear recognition of the value of computer programming for mathematics education. We use this theory and practice to suggest good techniques for and to raise questions about the use of AI in Mathematics Education.

A Bohemian matrix family is a set of matrices all of whose entries are drawn from a fixed, usually discrete and hence bounded, subset of a field of characteristic zero. Originally these were integers -- hence the name, from the acronym BOunded HEight Matrix of Integers (BOHEMI) -- but other kinds of entries are also interesting. Some kinds of questions about Bohemian matrices can be answered by numerical computation, but sometimes exact computation is better. In this paper we explore some Bohemian families (symmetric, upper Hessenberg, or Toeplitz) computationally, and answer some open questions posed about the distributions of eigenvalue densities.

In this paper, for certain type of structured { 0 , 1 , − 1 }–matrices, we give a complete description of the inner Bohemian inverses over any population containing the set { 0 , 1 , − 1 }. In addition, when the population is exactly { 0 , 1 , − 1 }, we provide explicit formulas for the number of inner Bohemian inverses of these type of matrices.