Search

Chebyshev polynomial spectral methods are very accurate, but are plagued by the cost and ill-conditioning of dense discretization matrices. Modified schemes, collectively known as “integration sparsification”, have mollified these problems by discretizing the highest derivative as a diagonal matrix. Here, we examine five case studies where the highest derivative diagonalization fails. Nevertheless, we show that Galerkin discretizations do yield banded matrices that retain most of the advantages of “integration sparsification”. Symbolic computer algebra greatly extends the reach of spectral methods. When spectral methods are implemented using exact rational arithmetic, as is possible for small truncation N in Maple, Mathematica and their ilk, roundoff error is irrelevant, and sparsification failure is not worrisome. When the discretization contains a parameter L , symbolic algebra spectral methods return, as answer to an eigenproblem, not discrete numbers but rather a plane algebraic curve defined as the zero set of a bivariate polynomial P ( λ , L ) ; the optimal approximations to the eigenvalues λ j are in the middle of the straight portions of the zero contours of P ( λ ; L ) where the isolines are parallel to the L axis.

When members of a classroom community routinely listen to one another and build on one another's ideas, not only do students learn and improve their cognitive and communication skills, but teacher and students develop a disposition to listen, think and talk together. However, such dialogic, classroom talk is rare. In this study we show how a teacher's epistemological commitment (that student ideas matter) combined with oracy practices (safe space for student talk; student ideas drive classroom talk; support multiple perspectives) realized through speculation and reasoning (S&R) words foster dialogic talk. We examined S&R words (think, would, might/maybe, if, so, but, how, why) in 1299 turns of talk in two lessons in one classroom of six 4–5th grade English Language Learners. Statistical discourse analysis showed that S&R words occurred more often during what we refer to as connect episodes (students made personal connections to the content), not after particular types of turns. Close discourse analysis showed how S&R word use cultivated a language of possibility and how management of classroom discourse modality promoted dialogic talk.

We analytically and numerically solve the PCP equation, u x x u y y = 1 , with homogeneous Dirichlet boundary conditions on the unit square. Chebyshev and Fourier spectral methods with low degree truncations yield moderate accuracy but the usual exponential rate of convergence of spectral methods is destroyed by the boundary singularities of the solution. In the sequel to this work, we will apply a variety of strategies including a change-of-coordinates and singular basis functions to recover spectral accuracy in spite of the boundary singularities. As preparation for this numerical study, we find explicit solutions to related problems to the two-dimensional PCP equation in a domain with a boundary that is an ellipse and the three-dimensional PCP equation in a cubic domain. We also analyze the boundary behavior of these solutions: all have complicated singularities with unbounded first derivatives.

We revisit the spectral solution of the Thomas–Fermi problem for neutral atoms, u r r − ( 1 / r ) u 3 / 2 = 0 on r ∈ [ 0 , ∞ ] with u ( 0 ) = 1 and u ( ∞ ) = 0 to illustrate some themes in solving differential equations when there are complications, and also to make improvements in our earlier treatment. By “complications” we mean features of the problem that either destroy the exponential accuracy of a standard Chebyshev series, or render the classic Chebyshev approach inapplicable. The Thomas–Fermi problem has four complications: (i) a semi-infinite domain r ∈ [ 0 , ∞ ] (ii) a square root singularity in u ( r ) at the origin (iii) a fractional power nonlinearity and (iv) asymptotic decay as r → ∞ that includes negative powers of r with fractional exponents. Our earlier treatment determined the slope at the origin to twenty-five decimal places, but no fewer than 600 basis functions were required to approximate a univariate solution that is everywhere monotonic, and all of the earlier tricks failed to recover an exponential rate of convergence in the truncation of the spectral series N, but only a high order convergence in negative powers of N. Here, using the coordinate z ≡ r to neutralize the square root singularity as before, we show that accuracy and rate of convergence are significantly improved by solving for the original unknown u ( r ) instead of the modified unknown v ( r ) = u ( r ) used previously. Without a further change of coordinate, a rational Chebyshev basis T L n ( z ; L ) yields twelve decimal place accuracy for the slope at the origin, u r ( 0 ) , with 70 basis functions and twenty-four decimal places with a truncation N = 100 . True exponential accuracy can be restored by using an appropriate change of coordinate, z = G ( Z ) , where G is some species of exponential. However, the various Chebyshev and Fourier series for the Thomas–Fermi function have “plural asymptotics”, that is, a n ∼ a intermediate ( n ) for 1 ≪ n ≪ n I but a n ∼ a far ( n ) for n ≫ n I for some positive constant n I . The “far-asymptotics” as n → ∞ are often of no practical significance. For this problem, the TL-with-sinh method, theoretically the best for huge n, is bedeviled with numerical ill-conditioning and its asymptotic superiority is realized only when the goal is at least forty decimal places of accuracy. This is absurd for engineering, but useful perhaps for benchmarking. To sixty decimal places, u r ( 0 ) = − 1.588071022611375312718684509423950109452746621674825616765677 . To nineteen digits after the decimal point, the constant in the Coulson–March asymptotic series is improved to F = 13.2709738480269351535 .

Purpose This study reviews contemporary management and follow-up of pediatric ovarian torsion. Methods This is a retrospective series of patients from birth to 19 years undergoing operative management of ovarian torsion from 2012 to 2016. Results We studied 43 girls who underwent 51 operations for ovarian torsion. The median age was 8.3 years. Ultrasound was utilized for diagnosis in 24/29 patients (83%) evaluated in a children's hospital. In contrast, computed tomography was used initially in 7 cases (50%) in children imaged at non–children's hospitals before transfer. Initial operation for ovarian torsion was completed laparoscopically in 38 (88%). Overall, ovarian preservation was performed in 37 (86%) patients, while 6 (13%) underwent oophorectomy. Indications for oophorectomy included 5 infants with in utero torsion and an 18-year-old with a suspected malignancy. In girls with acute ovarian torsion, the oophorectomy rate was reduced to 2%. Postoperatively, 1 patient developed a small bowel obstruction requiring operation after laparoscopic ovarian detorsion. Recurrent torsion occurred in 3 patients (7%). In total, 34 patients underwent postoperative ovarian imaging. A total of 25 (74%) had follicles visualized in the previously torsed ovary. Conclusion Ovarian-sparing operations for acute torsion are safe and result in ovarian salvage and preservation of follicular development in more than 70% of children and adolescents.

There is great interest in reaching a consensus method for in vitro digestion of foods. The COST INFOGEST network of researchers has developed a standardized protocol, which is appropriate for basic digestibility studies, but may be insufficient for more complex studies, such as micronutrient bioaccessibility. In the present study, the bioaccessibility of phytosterols (PS) was evaluated using the standardized and modified versions of the in vitro digestion protocol. The INFOGEST protocol was used without modification for 2% solutions of either free phytosterols (FPS) or esterified phytosterols (EPS) in soybean oil. The effect of two variables which have yet to be explored using the INFOGEST protocol was subsequently evaluated. Those variables are: (1) changing the mode of agitation from orbital shaking to head over heels (HOH) rotation and (2) adding pancreatic cholesterol esterase (CE) during the intestinal phase. Bioaccessibility of FPS was greater than that of EPS, and HOH rotation caused a further increase in bioaccessibility of FPS but not EPS. Addition of CE caused a slight, but significant improvement in bioaccessibility of EPS, though hydrolysis was incomplete. Our findings highlight a need for further evaluation of the standardized in vitro digestion protocol, especially for complex studies such as lipid bioaccessibility.

The nonlinear two-point boundary value problem (TPBVP for short) u x x + u 3 = 0 , u ( 0 ) = u ( 1 ) = 0 , offers several insights into spectral methods. First, it has been proved a priori that ? 0 1 u ( x ) d x = p / 2 . By building this constraint into the spectral approximation, the accuracy of N + 1 degrees of freedom is achieved from the work of solving a system with only N degrees of freedom. When N is small, generic polynomial system solvers, such as those in the computer algebra system Maple, can find all roots of the polynomial system, such as a spectral discretization of the TPBVP. Our second point is that floating point arithmetic in lieu of exact arithmetic can double the largest practical value of N . (Rational numbers with a huge number of digits are avoided, and eliminating M symbols like 2 and p reduces N + M -variate polynomials to polynomials in just the N unknowns.) Third, a disadvantage of an “all roots” approach is that the polynomial solver generates many roots – ( 3 N - 1 ) for our example – which are genuine solutions to the N -term discretization but spurious in the sense that they are not close to the spectral coefficients of a true solution to the TPBVP. We show here that a good tool for “root-exclusion” is calculating ? = ? n = 1 N b n 2 ; spurious roots have ? larger than that for the physical solution by at least an order of magnitude. The ? -criterion is suggestive rather than infallible, but root exclusion is very hard, and the best approach is to apply multiple tools with complementary failings.

The boundary value problem uuxx=2 appears in Cranes theory of laminar convection from a point source. We show that the solution is real only when |x|/2. On this interval, denoting the constants of integration by A and s, the general solution is AV([xs]/A) where the Crane function V is the parameter-free function V=exp({erfinv([2/])x}2) and erfinv(z) is the inverse of the error function. V(x) is weakly singular at both endpoints; its Chebyshev polynomial coefficients an decrease proportionally to 1/n3. Exponential convergence can be restored by writing V(x)=n=0a2nT2n(z[x]) where the mapping is z=arctanh(x/)L2+(arctanh(x/))2,=/2. Another option is singular basis functions. V(1x2/2){10.216log(1x2/2)} has a maximum pointwise error that is less 1/2000 of the maximum of the Crane function.

The Petrov-Galerkin ultraspherical polynomial/Chebyshev polynomial discretization of the highest derivative of a differential equation is a diagonal matrix. The same is true for Fourier-Galerkin discretizations. Nevertheless, the spectral discretizations of simple problems like u x x + q ( x ) u = f ( x ) are usually dense matrices. The villain is the "multiplication matrix", the Galerkin representation of a term like q ( x ) u ( x ) ; unfortunately, this part of the Galerkin matrix is dense. However, if the ODE coefficient q ( x ) has a Chebyshev or Fourier series that converges much more rapidly than u ( x ) , then it is possible to realize great cost savings at no loss of accuracy by truncating the full N ? N Galerkin matrix to a banded matrix where the bandwidth m ? N . One of our themes is that when the spectral series for q ( x ) and u ( x ) have similar rates of convergence, as is almost universal when a nonlinear equation is linearized for a Newton-Krylov iteration, such "accuracy lossless" truncation is impossible. Nonlinearity is but one of many causes of this sort of solution/coefficient "equiconvergence". When bandwidth truncation is possible, though, our second theme is to show that a modest amount of floating point operations and memory can be saved by an unsymmetric truncation in which the number of elements retained to the left of the main diagonal is roughly double the number kept to the right. Our second improvement is to replace the M -term spectral series for q ( x ) by its ( M / 2 ) / ( M / 2 ) Chebyshev-Pade rational approximation. This sometimes allow one to halve the matrix bandwidth, reducing the linear algebra costs by a factor of four.

Hitherto, as a tool for tracing all branches of nonlinear differential equations, resolution-increasing homotopy methods have been applied only to finite difference discretizations. However, spectral Galerkin algorithms typically match the error of fourth order differences with one-half to one-fifth the number of degrees of freedom N in one dimension, and a factor of eight to a hundred and twenty-five in three dimensions. Let u → N be the vector of spectral coefficients and R → N the vector of N Galerkin constraints. A common two-part procedure is to first find all roots of R → N ( u → N ) = 0 → using resultants, Groebner basis methods or block matrix companion matrices. (These methods are slow and ill-conditioned, practical only for small N .) The second part is to then apply resolution-increasing continuation. Because the number of solutions is an exponential function of N , spectral methods are exponentially superior to finite differences in this context. Unfortunately, u → N is all too often outside the domain of convergence of Newton’s iteration when N is increased to ( N + 1 ) . We show that a good option is the artificial parameter homotopy H → ( u → ; τ ) ≡ R → N + 1 ( u → ) − ( 1 − τ ) R → N + 1 ( u → N ) , τ ∈ [ 0 , 1 ] . Marching in small steps in τ , we proceed smoothly from the N -term to the N + 1 -term approximations.

The Charney problem, a second order ordinary differential equation eigenproblem on z ? 0 , ∞ with complex eigenvalues, is of great historical importance in meteorology and oceanography. Here, it is used as a testbed for several extensions of spectral methods. The first is to parameterize a plane curve which is singular at an endpoint, as very common in applications. The second stretch is to extend the Chebyshev tau method to compute eigenfunctions of the form M ( z ) + log ( z ) V ( z ) where M ( z ) and V ( z ) are entire functions and where the approximation interval is a line segment in the complex plane. Third, we offer a special procedure for finding the roots of a function which is not a polynomial, but rather the combination of a polynomial plus a logarithm multiplied by a second polynomial. Lastly, to resolve the very thin boundary layer of the regularized Charney problem, we combine a rational Chebyshev ( T L n ) pseudospectral method with a change of coordinate which is quadratic at the ground. Remarkably, best results are obtained by applying four boundary conditions even though the Charney problem is a differential equation of only second order.

Data are from a two-year ethnographic study of a diverse, urban classroom community of seven-and-eight-year-olds in the US. Research interest is in student uptake in their writing as a marker of learning. We want to know whether, and how, students who have routinely experienced exploratory talk about instructional points, take up these points in their writing a few weeks later. We examine a class set of writing for evidence. By student uptake we mean the appearance of a teaching point in individually written student stories. We apply nine uptake codes and report descriptive statistics for student stories (N=17) and take a qualitative look at 6 stories. We found widespread and consistent uptake of instructional points.

Fourier-Galerkin discretizations of PDEs can often be preconditioned by a matrix that is the union of a dense M × M matrix plus a large matrix that is either diagonal or banded with small bandwidth. This was suggested forty years ago by L.M. Delves. For our two-dimensional example, the preconditioned iteration converges geometrically fast to machine precision in less than ten iterations. If the residual of the partial differential equation is evaluated by Fast Fourier Transform (FFT), the cost scales almost linearly in the number of unknowns. On a 6144 × 6144 grid, the computation needed less than an hour in Matlab on a laptop. The dense block requires a total degree ordering; the FFT evaluation of the residual requires speedometer ordering of the Fourier coefficient matrix. We show that this dual bookkeeping is essential but not difficult.

The complete elliptic integral of the first kind arises in many applications. This article furnishes four different ways to compute the inverse of the elliptic integral. One motive for this study is simply that the author needed to compute the inverse integral for an application. Another is to develop a case study comparing different options for solving transcendental equations like those in the author’s book (Boyd, 2014). A third motive is to develop analytical approximations, more useful to theorists than mere numbers. A fourth motive is to provide robust “black box” software for computing this function. The first solution strategy is “polynomialization” which replaces the elliptic integral by an exponentially convergent series of Chebyshev polynomials. The transcendental equation becomes a polynomial equation which is easily solved by finding the eigenvalues of the Chebyshev companion matrix. (The numerically ill-conditioned step of converting from the Chebyshev to monomial basis is never necessary). The second approximation is a regular perturbation series, accurate where the modulus is small. The third is a power-and-exponential series that converges over the entire range parameter range, albeit only sub-exponentially in the limit of zero modulus. Lastly, Newton’s iteration is promoted from a local iteration to a global method by a Never-Failing Newton’s Iteration (NFNI) in the form of the exponential of the ratio of a linear function divided by another linear polynomial. A short Matlab implementation is provided, easily translatable into other languages. The Matlab/Newton code is recommended for numerical purposes. The other methods are presented because (i) all are broadly applicable strategies useful for other rootfinding and inversion problems (ii) series and substitutions are often much more useful to theorists than numerical software and (iii) the Never-Failing Newton’s Iteration was discovered only after a great deal of messing about with power series, inverse power series and so on.

Single-domain spectral methods have been largely restricted to tensor product bases on a tensor product grid. To break the "tensor barrier", we studied approximation in two idealized families of domains. One family is bounded by a "squircle", the zero isoline of B ( x , y ) = x 2 ? + y 2 ? - 1 . The boundary varies smoothly from a circle ? = 1 to the square ? = ∞ . The other family is bounded by a "perturbed quadrifolium", the plane algebraic curve ? ( x 2 + y 2 ) - ( ( x 2 + y 2 ) 3 - ( x 2 - y 2 ) 2 ) ; this varies smoothly from the singular, self-intersecting curve known as the quadrifolium to a circle as ? varies from zero to infinity. We compared two different bivariate polynomial bases, truncating by total degree. Zernike polynomials are natural for the disk; tensor products of Chebyshev polynomials are equally sensible for the square. Both yield an exponential rate of convergence for our non-tensor, neither disk-nor-square domains; indeed, the Chebyshev basis worked well for the disk and the Zernike polynomials were good for the square. The expected differences due to numerical ill-conditioning did not emerge, much to our surprise. The price for the nontensor domain was that hyperinterpolation was necessary, that is, least squares fitting with more interpolation constraints than unknowns. Denoting the number of interpolation points by P and the basis size by N, a ratio of P/N around two to three was optimum while P near one was very inacccurate. A uniform grid, truncated to include only those points within the squircle or other boundary curve, was satisfactory even without interpolation points on the boundary (although boundary points are a cost-effective improvement). Interpolation costs were greatly reduced by exploiting the invariance of the squircle-bounded and perturbed-quadrifolium domains to the eight element D4 dihedral group.

On an infinite grid with uniform spacing h, the cardinal basis Cj(x; h) for many spectral methods consists of translates of a "master cardinal function", Cj(x; h) = C(x/h - j). The cardinal basis satisfies the usual Lagrange cardinal condition, Cj(mh) = ?jm where ?jm is the Kronecker delta function. All such "shift-invariant subspace" master cardinal functions are of "localized-sinc" form in the sense that C(X) = sinc(X)s(X) for a localizer function s which is smooth and analytic on the entire real axis and the Whittaker cardinal function is sinc(X) ? sin?(πX)/(πX). The localized-sinc approximation to a general f(x) is f localized - sinc ( x ; h ) ? ? j = - ∞ ∞ f ( j h ) s ( x - j h ] / h ) sinc ( x - jh ] / h ) . In contrast to most radial basis function applications, matrix factorization is unnecessary. We prove a general theorem for the Fourier transform of the interpolation error for localized-sinc bases. For exponentially-convergent radial basis functions (RBFs) (Gaussians, inverse multiquadrics, etc.) and the basis functions of the Discrete Singular Convolution (DSC), the localizer function is known exactly or approximately. This allows us to perform additional error analysis for these bases. We show that the error is similar to that for sinc bases except that the localizer acts like a diffusion in Fourier space, smoothing the sinc error.

Bipolar coordinates provide an efficient cartography for a variety of geometries: the exterior of two disks or cylinders, a half-plane containing a disk, an eccentric annulus with a small disk offset from the center of an outer boundary that is a large circle, and so on. A pseudospectral method that employs a tensor product basis of Fourier functions in the cyclic coordinate ? and Chebyshev polynomials in the quasi-radial coordinate ? gives easy-to-program spectral accuracy. We show, however, that as the inner disk becomes more and more offset from the center of the outer boundary circle, the grid is increasingly non-uniform, and the rate of exponential convergence increasingly slow. One-dimensional coordinate mappings significantly reduce the non-uniformity. In spite of this non-uniformity, the Chebyshev-Fourier method is quite effective in an idealized model of the wind-driven ocean circulation, resolving both internal and boundary layers. Bipolar coordinates are also a good starting point for solving problems in a domain which is not one of the "bipolar-compatible" domains listed above, but is a sufficiently small perturbation of such. This is illustrated by applying boundary collocation with bipolar harmonics to solve Laplace's equation in a perturbed eccentric annulus in which the disk-shaped island has been replaced by an island bounded by an ellipse. Similarly a perturbed bipolar domain can be mapped to an eccentric annulus by a smooth change of coordinates.

Vortex blob methods approximate a flow as a sum of many small vortices of Gaussian shape, and adaptively move the vortex centers with the current. Gaussian radial basis functions (RBFs) do exactly the same. However, RBFs solve an exact interpolation problem - expensive but accurate - while vortex methods sacrifice accuracy through quasi-interpolation for the absence of a matrix inversion. We show that vortex-RBF algorithms with spectral accuracy are stable for flows on the sphere. The version in Eulerian coordinates is fast; the fully-Lagrangian variant is much slower for a given basis size N, but is highly adaptive for advection-dominated flows. Both versions are excellent for small-to-medium-N problems - N up to 10,000, say, where N is the number of RBF grid points/vortex blobs. Neither is good for large N problems because the cost of the Eulerian model scales as N 2 per timestep while the Lagrangian vortex-RBF method scales as N 3 . The slow-Lagrangian scheme is unique among vortex methods in being genuinely (like its Eulerian sibling) a spectrally-accurate method: for laminar flows the error falls exponentially fast with N.

The electrode in a cochlear implant (CI) system is a key factor in hearing performance as it is the interface between the device and the auditory pathway of the recipient. The first CI electrodes were straight and thus adopted a lateral wall position. Subsequent innovations include: perimodiolar electrodes designed to lie adjacent to the modiolar wall and thus to provide more spatially-focused stimulation of the spiral ganglion cells; shorter atraumatic straight electrodes for combined electric and acoustic (hybrid) stimulation. This paper explores the relative merits of straight and perimodiolar electrodes in the search for the optimal electrode design with reference to electrodes from Cochlear(®).

Asymptotic approximations show that the lowest modes of the prolate spheroidal wave equation are concentrated with an O ( 1 / c ) length scale where c is the "bandwidth" parameter of the prolate differential equation. Accurate computation of the ground state eigenvalue by the long-known Legendre-Galerkin method requires roughly 3.8 c Legendre polynomials. Some studies have therefore applied a grid with 20,000 points in conjunction with high order finite differences to achieve c = 10 7 . Here, we show that by adaptively applying either Hermite functions or Gaussian radial basis functions (RBFs), it is never necessary to use more than eighty degrees of freedom to calculate the lowest dozen eigenvalues of each symmetry class. For small c, the eigenmodes are not confined to a small portion of the domain ? ? - π / 2 , π / 2 in latitude, but are global. We show that by periodizing the basis functions via imbricate series, it is possible to apply Hermite and RBFs even in the limit c ? 0 . (The Legendre method is probably a little more efficient in this limit since the prolate functions tend to Legendre polynomials in this limit.) A "sideband truncation" restricts the discretization to a small block taken from the large Hermite-Galerkin matrix. We show that sideband truncation with blocks as small as 5 × 5 is a very efficient way to compute high order modes. In an appendix, we prove a rigorous convergence theorem for the periodized Hermite expansion.

We extend the theory of periodized RBFs. We show that the imbricate series that define the Periodic Gaussian (PGA) and Sech (PSech) basis functions are Jacobian theta functions and elliptic functions “dn”, respectively. The naive periodization fails for the Multiquadric and Inverse Multiquadric RBFs, but we are able to define periodic generalizations of these, too, by proving and exploiting a generalization of the Poisson Summation Theorem. Although applications of periodic RBFs are mostly left for another day, we do illustrate the flaws and potential by solving the Mathieu eigenproblem on both uniform and highly-adapted grids. The terms of a Fourier basis can be grouped into four classes, depending upon parity with respect to both the origin and x = π / 2 , and so, too, the Mathieu eigenfunctions. We show how to construct symmetrized periodic RBFs and illustrate these by solving the Mathieu problem using only the periodic RBFs of the same symmetry class as the targeted eigenfunctions. We also discuss the relationship between periodic RBFs and trigonometric polynomials with the aid of an explicit formula for the nonpolynomial part of the Periodic Inverse Quadratic (PIQ) basis functions. We prove that the rate of convergence for periodic RBFs is geometric, that is, the error can be bounded by exp ⁡ ( − N μ ) for some positive constant μ. Lastly, we prove a new theorem that gives the periodic RBF interpolation error in Fourier coefficient space. This is applied to the “spectral-plus” question. We find that periodic RBFs are indeed sometimes orders of magnitude more accurate than trigonometric interpolation even though it has long been known that RBFs (periodic or not) reduce to the corresponding classical spectral method as the RBF shape parameter goes to 0. However, periodic RBFs are “spectral-plus” only when the shape parameter α is adaptively tuned to the particular f ( x ) being approximated and even then, only when f ( x ) satisfies a symmetry condition.

The Fourier Transform of a quartic Gaussian, @F(k)=@!"-"~^~exp(ikx)exp(-x^4)dx, is important in the theory of radial basis functions with @f(r)=exp(-r^4). We show that the transform is the exact sum of two generalized hypergeometric functions "0F"2. This implies that @F(k) is an entire function and also gives an analytical form for the power series coefficients of arbitrary degree. Using the method of steepest descent for integrals, we derive the first five terms in an asymptotic series in powers of k^-^4^/^3. Through a simple hyperasymptotic analysis, we show that the magnitude of the error of the optimally-truncated series [''superasymptotic'' error] is Oexp-0.818k^4^/^3. The lowest order approximation, @F(k)~2^7^/^6@p/3k^-^1^/^3exp(-(2^1^/^33/16)k^4^/^3)cos((3^3^/^22^1^/^3/16)k^4^/^3-@p/6) shows that the transform @F(k) is an exponentially-decaying oscillation. We show that the steepest descent method [2,27] (Bender and Orszag, 1978; Miller, 2006) [26,33] (Lauwerier, 1966; Ursell, 1965) can be applied to the transform of exp(-Ax^2^n) for general integer n and give the first two terms. Another theme is that modern computer technology has simplified the steepest descent method. The steepest descent paths can be plotted by merely displaying the contours of the imaginary part of the phase function @J(t) where z@J(t) is the logarithm of the integrand, and z is the large parameter (here proportional to k^4^/^3). The task of inverting w(t) to obtain a series in powers of w for the metric factor dt/dw after the usual change of integration variable can now be done to high order by a few lines of an algebraic manipulation language as illustrated by a short table with the complete Maple code.

We revisit early models of steady western boundary currents [Gulf Stream, Kuroshio, etc.] to explore the role of irregular coastlines on jets, both to advance the research frontier and to illuminate for education. In the framework of a steady-state, quasigeostrophic model with viscosity, bottom friction and nonlinearity, we prove that rotating a straight coastline, initially parallel to the meridians, significantly thickens the western boundary layer. We analyze an infinitely long, straight channel with arbitrary orientation and bottom friction using an exact solution and singular perturbation theory, and show that the model, though simpler than Stommel's, nevertheless captures both the western boundary jet (“Gulf Stream”) and the “orientation effect”. In the rest of the article, we restrict attention to the Stommel flow (that is, linear and inviscid except for bottom friction) and apply matched asymptotic expansions, radial basis function, Fourier–Chebyshev and Chebyshev–Chebyshev pseudospectral methods to explore the effects of coastal geometry in a variety of non-rectangular domains bounded by a circle, parabolas and squircles. Although our oceans are unabashedly idealized, the narrow spikes, broad jets and stationary points vividly illustrate the power and complexity of coastal control of western boundary layers.

A domain is invariant under the eight-element D 4 symmetry group if it is unchanged by reflection with respect to the x and y axes and also the diagonal line x = y . Previous treatments of group theory for spectral methods have generally demanded a semester?s worth of group theory. We show this is unnecessary by providing explicit recipes for creating grids, etc. We show how to decompose an arbitrary function into six symmetry-invariant components, and thereby split the interpolation problem into six independent subproblems. We also show how to make symmetry-invariant basis functions from products of Chebyshev polynomials, from Zernike polynomials and from radial basis functions (RBFs) of any species. These recipes are completely general, and apply to any domain that is invariant under the dihedral group D 4 . These concepts are illustrated by RBF pseudospectral solutions of the Poisson equation in a domain bounded by a squircle, the square-with-rounded corners defined by x 2 ? + y 2 ? - 1 = 0 where here ? = 2 . We also apply Chebyshev polynomials to compute eigenmodes of the Helmholtz equation on the square and show each mode belongs to one and only one of the six D 4 classes.

A frequent claim in the radial basis function (RBF) literature is that with a wise choice of the ''shape parameter'' @a, spectrally accurate RBFs can be ''spectral-plus'' in the sense that the RBF approximation is orders of magnitude more accurate than the corresponding spectral method. This is surprising because in the so-called ''flat limit'', @a->0, RBF methods provably reduce to a standard spectral method. On an infinite uniform grid, for example, Gaussian, hyperbolic secant, inverse quadratic, and multiquadric and inverse multiquadric RBFs reduce to the classical sinc basis. Using the error theorem proved in Part 1, we show that indeed it is possible for RBFs to be greatly superior to the sinc basis. However, rather special conditions are required for this superiority. The generic case is that RBFs are less accurate than the sinc method for nonzero values of the shape parameter.

Analysis of the Sea Surface Height (SSH) from satellite altimeters has shown that equatorially trapped Rossby waves exhibit asymmetric cross-equatorial structures; their northern extrema are much larger in magnitude than their southern counterparts. Such asymmetry is inconsistent with the classical theory for the first baroclinic, first meridional equatorially trapped Rossby mode, which predicts that SSH and zonal velocity are symmetric in latitude and the meridional velocity is latitudinally antisymmetric ( Matsuno, 1966 ). Chelton et al. (2003) attributed the observed asymmetry to the mean-shear-induced modifications of first meridional mode Rossby waves. The present paper examines nonlinear rectification of cross-equatorial wave structures in the presence of different zonal mean currents. Nonlinear traveling Rossby waves embedded in shears are calculated numerically in a 1.5-layer model. Nonlinearity is shown to increase the cross-equatorial asymmetry substantially making the northern extrema even more pronounced. However, nonlinearity only slightly increases the magnitude of the westward phase speed.

A trigonometric polynomial is a truncated Fourier series of the form f N ( t ) ≡ ∑ j = 0 N a j cos ( jt ) + ∑ j = 1 N b j sin ( jt ) . It has been previously shown by the author that zeros of such a polynomial can be computed as the eigenvalues of a companion matrix with elements which are complex valued combinations of the Fourier coefficients, the “CCM” method. However, previous work provided no examples, so one goal of this new work is to experimentally test the CCM method. A second goal is introduce a new alternative, the elimination/Chebyshev algorithm, and experimentally compare it with the CCM scheme. The elimination/Chebyshev matrix (ECM) algorithm yields a companion matrix with real-valued elements, albeit at the price of usefulness only for real roots. The new elimination scheme first converts the trigonometric rootfinding problem to a pair of polynomial equations in the variables ( c , s ) where c ≡ cos ( t ) and s ≡ sin ( t ) . The elimination method next reduces the system to a single univariate polynomial P ( c ) . We show that this same polynomial is the resultant of the system and is also a generator of the Groebner basis with lexicographic ordering for the system. Both methods give very high numerical accuracy for real-valued roots, typically at least 11 decimal places in Matlab/IEEE 754 16 digit floating point arithmetic. The CCM algorithm is typically one or two decimal places more accurate, though these differences disappear if the roots are “Newton-polished” by a single Newton’s iteration. The complex-valued matrix is accurate for complex-valued roots, too, though accuracy decreases with the magnitude of the imaginary part of the root. The cost of both methods scales as O ( N 3 ) floating point operations. In spite of intimate connections of the elimination/Chebyshev scheme to two well-established technologies for solving systems of equations, resultants and Groebner bases, and the advantages of using only real-valued arithmetic to obtain a companion matrix with real-valued elements, the ECM algorithm is noticeably inferior to the complex-valued companion matrix in simplicity, ease of programming, and accuracy.

We solve the Thomas-Fermi problem for neutral atoms, u"y"y-(1/y)u^3^/^2=0 on y@?[0,~] with u(0)=1 and u(~)=0, using rational Chebyshev functions TL"n(y;L) to illustrate some themes in solving differential equations on a semi-infinite interval. L is a user-choosable numerical parameter. The Thomas-Fermi equation is singular at the origin, giving a TL convergence rate of only fourth order, but this can be removed by the change of variables, z=y with v(z)=u(y(z)). The function v(z) decays as z->~ with a term in z^-^3, which is consistent with a geometric rate of convergence. However, the asymptotic series has additional terms with irrational fractional powers beginning with z^-^4^.^5^4^4. In spite of the faster spatial decay, the irrational powers degrade the convergence rate to slightly larger than tenth order. This vividly illustrates the subtle connection between the spatial decay of u(x) and the decay-with-degree of its rational Chebyshev series. The TL coefficients a"n(L) are hostages to a tug-of-war between a singularity on the negative real axis, which gives a geometric rate of convergence that slows with increasingL, and the slow inverse power decay for large z, which gives quasi-tenth order convergence with a proportionality constant that decreasesinversely as a power of L. For L=2, we can approximate u"y(0) (=v"z"z(0)) to 1 part in a million with a truncation N of only 20. L=64 and N=600 gives u"y(0)=-1.5880710226113753127186845, correct to 25 decimal places.

We catalogue the numerical properties of approximations using two novel types of radial basis functions @f(r). The QG species is a basis of exponentials of quartic argument: f(x)~f^R^B^F(x;@a,h)=@?"j"="1^Na"jexp(-[@a/h]^4(x-x"j)^4) where the x"j are the RBF centers and also the interpolation points. We show that Quartic Gaussian RBFs fail at many discrete values of the shape parameter @a. We show through a detailed analysis that these singularities are directly related to zeros of Q(k), the Fourier Transform of exp(-x^4). If we reverse the roles and take Q(x) as the RBF, all difficulties disappear because these IQG RBFs have a Fourier transform which is nonnegative for all real k. We explain that although the Quartic-Gaussian exp(-x^4) is positive definite in the physics/dynamical systems sense of being zero-free and nonnegative, it lacks the crucial property of being positive definition in the RBF/analysis sense.

A function is symmetric with respect to a point x = L if f ( x + L ) = f ( - x + L ) for all x and similarly is antisymmetric if f ( x + L ) = - f ( - x + L ) . A function which is either symmetric or antisymmetric is said to be of “definite parity” with respect to L . The sines and cosines of a Fourier series have definite parity with respect to two points; all cosines are symmetric with respect to the origin while all sines are antisymmetric with respect to x = 0 ; cos ( 2 nx ) and sin ( [ 2 n + 1 ] x ) for integral n are also symmetric with respect to x = π / 2 while all other Fourier basis functions are antisymmetric with respect to the same point. Such symmetries can be exploited in numerical calculations; for example, computing the angular Mathieu functions using N basis functions can be split into four uncoupled eigenproblems each of dimension N / 4 . It is natural to ask: Are there other classes of functions with similar symmetries? Using concepts from computer graphics, we prove that all functions which are symmetric with respect to two points separated by a distance L must be spatially periodic with period 4 L . We also prove that the only function which is of definite parity with respect to three distinct points must be a constant. These theorems define parity in the usual sense of a global property such that even parity with respect to the origin means f ( x ) = f ( - x ) for all x ∈ [ - ∞ , ∞ ] . We construct counterexamples to both theorems that are functions with local parity, that is, symmetry which applies only for a finite interval in x .

Fourier series often need to be generalized by appending a linear polynomial to the usual series of sines and cosines. The integral of a trigonometric polynomial is one example; another is a time series of climate data where the periodic oscillations of the diurnal and annual cycles are accompanied by a non-periodic trend (global warming). Stock market averages fluctuate about a generally upward trend. Such non-periodic variations with time are commonly called “secular trends”. We borrow “secular” to label a truncated Fourier series plus a linear trend as a “linear”, secular trigonometric polynomial. Standard Fourier rootfinding methods are wrecked by the extra, nonperiodic term. Therefore, we introduce a new algorithm for computing the zeros of a Fourier polynomial-with-secular-trend. First, we expand the linear secular trigonometric polynomial f N ( t ) as a truncated Chebyshev series. Because of the special structure, it is easy to calculate a problem-dependent truncation M such that the error of the truncated Chebyshev series is guaranteed to be less than a user-specified tolerance. We then find the roots of the truncated Chebyshev series as the eigenvalues of the Chebyshev companion matrix. This computes all roots, but we explain why the method is not reliable for complex-valued roots unless these are close to the real axis. No a priori information is required of the user except the coefficients of the linear secular trigonometric polynomial. Numerical examples show that 13 decimal place accuracy for real roots is typical.

The incomplete elliptic integral of the second kind, E ( sin ( T ) , m ) ≡ ∫ 0 T 1 - m sin 2 ( T ′ ) dT ′ where m ∈ [ 0 , 1 ] is the elliptic modulus, can be inverted with respect to angle T by solving the transcendental equation E ( sin ( T ) ; m ) - z = 0 . We show that Newton’s iteration, T n + 1 = T n - E ( sin ( T ) ; m ) - z 1 - m sin 2 ( T ) , always converges to T ( z ; m ) = E - 1 ( z ; m ) within a relative error of less than 10 −10 in three iterations or less from the first guess T 0 ( z , m ) = π / 2 + r ( θ - π / 2 ) where, defining ζ ≡ 1 - z / E ( 1 ; m ) , r = ( 1 - m ) 2 + ζ 2 and θ = atan ( ( 1 - m ) / ζ ) . We briefly discuss three alternative initialization strategies: “homotopy” initialization [ T 0 ( z , m ) ≡ ( 1 - m ) ( z ; 0 ) + mT ( z , m ) ( m ; 1 ) ], perturbation series (in powers of m ), and inversion of the Chebyshev interpolant of the incomplete elliptic integral. Although all work well, and are general strategies applicable to a very wide range of problems, none of these three alternatives is as efficient as the empirical initialization, which is completely problem-specific. This illustrates Tai Tsun Wu’s maxim “usefulness is often inversely proportional to generality”.

Computation of solitons of the cubically-nonlinear Benjamin-Ono equation is challenging. First, the equation contains the Hilbert transform, a nonlocal integral operator. Second, its solitary waves decay only as O(1/|x|^2). To solve the integro-differential equation for waves traveling at a phase speed c, we introduced the artificial homotopy H(u"X"X)-c u+(1-@d)u^2+@du^3=0, @d@?[0,1] and solved it in two ways. The first was continuation in the homotopy parameter @d, marching from the known Benjamin-Ono soliton for @d=0 to the cubically-nonlinear soliton at @d=1. The second strategy was to bypass continuation by numerically computing perturbation series in @d and forming Pade approximants to obtain a very accurate approximation at @d=1. To further minimize computations, we derived an elementary theorem to reduce the two-parameter soliton family to a parameter-free function, the soliton symmetric about the origin with unit phase speed. Solitons for higher order Benjamin-Ono equations are also computed and compared to their Korteweg-deVries counterparts. All computations applied the pseudospectral method with a basis of rational orthogonal functions invented by Christov, which are eigenfunctions of the Hilbert transform.

The Benjamin–Ono equation is especially challenging for numerical methods because (i) it contains the Hilbert transform, a nonlocal integral operator, and (ii) its solitary waves decay only as O(1/|x|2). We compare three different spectral methods for solving this one-space-dimensional equation. The Fourier pseudospectral method is very fast through use of the Fast Fourier Transform (FFT), but requires domain truncation: replacement of the infinite interval by a large but finite domain. Such truncation is unnecessary for a rational basis, but it is simple to evaluate the Hilbert Transform only when the usual rational Chebyshev functions TBn(x) are replaced by their cousins, the Christov functions; the FFT still applies. Radial basis functions (RBFs) are slow for a given number of grid points N because of the absence of a summation algorithm as fast as the FFT; because RBFs are meshless, however, very flexible grid adaptation is possible.

Through numerical experiments, we examine the condition numbers of the interpolation matrix for many species of radial basis functions (RBFs), mostly on uniform grids. For most RBF species that give infinite order accuracy when interpolating smooth f(x)-Gaussians, sech's and Inverse Quadratics-the condition number @k(@a,N) rapidly asymptotes to a limit @k"a"s"y"m"p(@a) that is independent of N and depends only on @a, the inverse width relative to the grid spacing. Multiquadrics are an exception in that the condition number for fixed @a grows as N^2. For all four, there is growth proportional to an exponential of 1/@a (1/@a^2 for Gaussians). For splines and thin-plate splines, which contain no width parameter, the condition numbers grows asymptotically as a power of N-a large power as the order of the RBF increases. Random grids typically increase the condition number (for fixed RBF width) by orders of magnitude. The quasi-random, low discrepancy Halton grid may, however, have a lower condition number than a uniform grid of the same size.

Approximating a smooth function from its values f(x"i) at a set of evenly spaced points x"i through P-point polynomial interpolation often fails because of divergence near the endpoints, the ''Runge Phenomenon''. This report shows how to achieve an error that decreases exponentially fast with P by means of polynomial interpolation on N"s subdomains where N"s increases with P. We rigorously prove that in the limit both N"s and M, the degree on each subdomain, increase simultaneously, the approximation error converges proportionally to exp(-constantPlog(P)). Thus, division into ever-shrinking, ever-more-numerous subdomains is guaranteed to defeat the Runge Phenomenon in infinite precision arithmetic. (Numerical ill-conditioning is also discussed, but is not a great difficulty in practice, though not insignificant in theory.) Although a Chebyshev grid on each subdomain is well known to be immune to the Runge Phenomenon, it is still interesting, and the same methodology can be applied as to a uniform grid. When a Chebyshev grid is used on each subdomain, there are two regimes. If c is the distance from the middle of the interval [-1,1] to the nearest singularity of f(x) in the complex plane, then when cN"s@?1, the error is proportional to exp(-cP), independent of the number of subdomains. When cN"s@?1, the rate of convergence slows to exp(-constantPlog(P)), the same as for equispaced interpolation. However, the Chebyshev multidomain error is always smaller than the equispaced multidomain error.

We compare seven different strategies for computing spectrally-accurate approximations or differential equation solutions in a disk. Separation of variables for the Laplace operator yields an analytic solution as a Fourier-Bessel series, but this usually converges at an algebraic (sub-spectral) rate. The cylindrical Robert functions converge geometrically but are horribly ill-conditioned. The Zernike and Logan-Shepp polynomials span the same space, that of Cartesian polynomials of a given total degree, but the former allows partial factorization whereas the latter basis facilitates an efficient algorithm for solving the Poisson equation. The Zernike polynomials were independently rediscovered several times as the product of one-sided Jacobi polynomials in radius with a Fourier series in ?. Generically, the Zernike basis requires only half as many degrees of freedom to represent a complicated function on the disk as does a Chebyshev-Fourier basis, but the latter has the great advantage of being summed and interpolated entirely by the Fast Fourier Transform instead of the slower matrix multiplication transforms needed in radius by the Zernike basis. Conformally mapping a square to the disk and employing a bivariate Chebyshev expansion on the square is spectrally accurate, but clustering of grid points near the four singularities of the mapping makes this method less efficient than the rest, meritorious only as a quick-and-dirty way to adapt a solver-for-the-square to the disk. Radial basis functions can match the best other spectral methods in accuracy, but require slow non-tensor interpolation and summation methods. There is no single "best" basis for the disk, but we have laid out the merits and flaws of each spectral option.

We discuss the numerical computation of the cosine lemniscate function and its inverse, the lemniscate integral. These were previously studied by Bernoulli, Euler, Gauss, Abel, Jacobi and Ramanujan. We review general elliptic formulas for this special case and provide some novelties. We show that a Fourier series by Ramanujan converges twice as fast as the standard elliptic cosine Fourier series specialized to this case. The so-called imbricate series, however, converges geometrically fast over the entire complex plane. We derive two new expansions. The rational-plus-Fourier series converges much faster than Ramanujan's: for real z: each term is asymptotically 12,400 times smaller than its immediate predecessor: coslem(z)=4B{q(1-q)cos(Bz)/[(1+q)^2-4qcos^2(Bz)]+@?"n"="1^~q^n^-^1^/^2{1/(1+q^2^n^-^1)-1}cos((2n-1)Bz)} where q=exp(-@p) is the elliptic nome, K~1.85... is the complete elliptic integral of the first kind for a modulus m=1/2 and B=@p/(K2). The rational imbricate series is uniformly valid over the complex plane, but converges twice as fast as the sech-imbricate series: coslem(z)=4Bq(1-q)@?"j"="-"~^~q(1-q)cos(Bz)/{(1+q)^2-4qcos^2(B[z-jPi])} where P=(4/2)K is the period in both the real and imaginary directions. We devise a new approximation for the lemniscate integral for real argument as the arccosine of a Chebyshev series and show that 17 terms yield about 15 digits of accuracy. For complex argument, we show that the lemniscate integral can be found to near machine precision (assumed as sixteen decimal digits) by computing the roots of a polynomial of degree thirteen. Alternatively, Newton's iteration converges in three iterations with an initialization, accurate to four decimal places, that is the chosen root of a cubic equation.

Radial basis function (RBF) interpolants have become popular in computer graphics, neural networks and for solving partial differential equations in many fields of science and engineering. In this article, we compare five different species of RBFs: Gaussians, hyperbolic secant (sech's), inverse quadratics, multiquadrics and inverse multiquadrics. We show that the corresponding cardinal functions for a uniform, unbounded grid are all approximated by the same function: C(X)~(1/(?))sin(πX)/sinh (πX/?) for some constant ?(α) which depends on the inverse width parameter ("shape parameter") α of the RBF and also on the RBF species. The error in this approximation is exponentially small in 1/α for sech's and inverse quadratics and exponentially small in 1/α2 for Gaussians; the error is proportional to α4 for multiquadrics and inverse multiquadrics. The error in all cases is small even for α~O(1).These results generalize to higher dimensions. The Gaussian RBF cardinal functions in any number of dimensions d are, without approximation, the tensor product of one dimensional Gaussian cardinal functions: C d ( x 1 , x 2 ? , x d ) = ? j = 1 d C ( x j ) . For other RBF species, we show that the two-dimensional cardinal functions are well approximated by the products of one-dimensional cardinal functions; again the error goes to zero as α?0. The near-identity of the cardinal functions implies that all five species of RBF interpolants are (almost) the same, despite the great differences in the RBF ?'s themselves.

We have applied the finite differences method to the study of a pair of isospectral heterogeneous domains, first introduced in Ref. [1]. We show that Richardson and Padé-Richardson extrapolations can be used (as in the homogeneous case) to obtain very precise approximations to the lowest eigenvalues. We have found that the first few exponents of the asymptotic series for the finite difference eigenvalues are unchanged with from the homogeneous case. Additionally, we have improved the previous best estimates for the case of homogeneous isospectral domains, obtaining 10 extra correct digits for the fundamental mode (and similar results for the other eigenvalues), with respect to the best result previously available. Keywords: Finite differences, Extrapolation, Isospectral domains

Artificial viscosity is a common device for stabilizing flows with shocks and fronts. The computational diffusion smears the frontal zone over a small distance μ where μ is chosen so that the discretization has a couple of grid points in the front, and thus is able to resolve the shock. Spectral element methods use a Legendre spectral viscosity whose effect is to damp the coefficient of P n ( x ) by some amount that depends only on the degree n of the Legendre polynomial. Legendre viscosity is better than ordinary diffusion because it does not require spurious boundary conditions, does not increase the temporal stiffness of the differential equations, and can be applied locally on an element-by-element basis. Unfortunately, Legendre diffusion is equivalent to a diffusion with a spatially-varying coefficient that goes to zero at the boundaries. Using the simplest example, one in which the second derivative of Burgers equation is replaced by the Legendre operator to give the “Legendre–Burgers” equation, u t + uu x = ν [(1 − x 2 ) u x ] x , we show that the width of the computational front can similarly tend to zero at the endpoints, causing a numerical catastrophe.

Background Cafe-au-lait macules (CALMs) in neurofibromatosis type 1 (NF1) are an early and accessible phenotype in NF1, but have not been extensively studied. Objective We sought to more fully characterize the phenotype of CALMs in patients with NF1. Methods In all, 24 patients with a diagnosis of NF1 confirmed through clinical diagnosis or molecular genetic testing were recruited from patients seen in the genetics department at the University of Alabama at Birmingham. CALM locations were mapped using standard digital photography. Pigment intensity was measured with a narrowband spectrophotometer, which estimates the relative amount of melanin based on its absorption of visible light. The major response was defined as the difference between the mean melanin from the CALM and the mean melanin from the surrounding skin. The major response for each spot was compared with spots within an individual and across individuals in the study population. Results There was significant variability of the major response, primarily attributable to intrapersonal variability (48.4%, P P Limitations The study was performed on a small population of patients and the method has not yet been used extensively for this purpose. Conclusions CALMs vary in pigment intensity not only across individuals, but also within individuals and this variability was unrelated to sun exposure. Further studies may help elucidate the molecular basis of this finding, leading to an increased understanding of the pathogenesis of CALMs in NF1.

Radial basis function (RBF) interpolation is a ''meshless'' strategy with great promise for adaptive approximation. One restriction is ''error saturation'' which occurs for many types of RBFs including Gaussian RBFs of the form @f(x;@a,h)=exp([email protected]^2(x/h)^2): in the limit h->0 for fixed @a, the error does not converge to zero, but rather to E"S(@a). Previous studies have theoretically determined the saturation error for Gaussian RBF on an infinite, uniform interval and for the same with a single point omitted. (The gap enormously increases E"S(@a).) We show experimentally that the saturation error on the unit interval, [email protected]?[-1,1], is about 0.06exp(-0.47/@a^2)@[email protected]?"~ - huge compared to the O([email protected]/@a^2)exp([email protected]^2/[[email protected]^2]) saturation error for a grid with one point omitted. We show that the reason the saturation is so large on a finite interval is that it is equivalent to an infinite grid which is uniform except for a gap of many points. The saturation error can be avoided by choosing @[email protected]?1, the ''flat limit'', but the condition number of the interpolation matrix explodes as O(exp(@p^2/[[email protected]^2])). The best strategy is to choose the largest @a which yields an acceptably small saturation error: If the user chooses an error tolerance @d, then @a"o"p"t"i"m"u"m(@d)=1/-2log(@d/0.06).

Radial basis functions are popular for interpolation on a scattered or irregular grid. However, theory for an irregular grid is mostly limited to proofs of convergence. Here, we present theory and numerical experiments for two specific cases. The first is an otherwise uniform grid of spacing h in which one point is shifted by an amount sh. The second is a uniform grid with one point omitted. We discuss Gaussian, hyperbolic secant, and inverse quadratic RBFs.

Gaussian radial basis functions (RBFs) on an infinite interval with uniform grid pacing h are defined by @f(x;@a,h)=exp(-[@a^2/h^2]x^2). The only significant numerical parameter is @a, the inverse width of the RBF functions relative to h. In the limit @a->0, we demonstrate that the coefficients of the interpolant of a typical function f(x) grow proportionally to exp(@p^2/[4@a^2]). However, we also show that the approximation to the constant f(x)=1 is a Jacobian theta function whose coefficients do not blow up as @a->0. The subtle interplay between the complex-plane singularities of f(x) (the function being approximated) and the RBF inverse width parameter @a are analyzed. For @a~1/2, the size of the RBF coefficients and the condition number of the interpolation matrix are both no larger than O(10^4) and the error saturation is smaller than machine epsilon, so this @a is the center of a ''safe operating range'' for Gaussian RBFs.

When diagonal values are missing or excluded, MINRES is a natural continuous model for the core/periphery structure of a symmetric social network matrix. Symmetric models, however, are not so useful when dealing with asymmetric data. Singular value decomposition (SVD) is a natural choice to model asymmetry, but this method also requires the presence of diagonal values. In this paper we offer an alternative, more general, approach to continuous core/periphery structures, the minimum residual singular value decomposition (MINRES/SVD), where each node in the network receives two indices, an “in-coreness” and an “out-coreness.” The algorithm for computing these coreness vectors is a least squares computation similar to, but distinct from the SVD, again because of the missing diagonal values. And in contrast to the standard, symmetric MINRES algorithm, we can more accurately model asymmetric matrices. This allows us to distinguish, for example, countries in the world economy that are more in the exporting core than they are in the importing core. We propose two nested PRE (proportional reduction of error) measures of fit: (1) the PRE from the MINRES vector with respect to the data and (2) the PRE of the product of the two MINRES/SVD vectors. Applying the resulting method to citations between journals and to international trade in clothing, we illustrate insights gained from being able to model asymmetrical flow patterns. Finally, two permutation tests are introduced to test independently for the MINRES and MINRES/SVD results.