User j. m. - MathOverflow

Answer by J. M. for Is it possible to generalize functions like $x^y, \ln x, \sin x, \arctan x$ to surreal numbers or surcomplex numbers?

2013-05-30T18:53:44Z

(Too long for a comment.)

It would seem to me that one other possibility for defining a transcendental function over the surcomplexes might be to use a suitable modification of the Cauchy integral formula; since you can reciprocate surcomplex numbers, one could then consider

$$f(z)=\frac1{2\pi i}\oint_\gamma f(t) (t-z)^{-1}\mathrm dt$$

where $\gamma$ is some suitable anticlockwise contour. I understand that this is on the surface a naïve proposal, and I would be interested in hearing how this might break down.

Answer by J. M. for Evaluating a limit similar to the Euler constant

2012-08-07T01:19:05Z

As already mentioned by Fedor and Igor, you can ignore the logarithmic term since it zeroes out at $\infty$, and you can just concentrate on the series

$$-\frac{i}{2}\sum_{k=1}^\infty \frac1{ik+k^{3/2}}$$

Using Laplace transform techniques, your sum can be transformed into the integral

$$-\frac12-\frac{i}{\sqrt\pi}\int_0^\infty \frac{F(\sqrt{u})}{\exp\,u-1}\mathrm du$$

where $F(z)$ is Dawson's integral.

I don't know of a closed form for this integral, but Mathematica easily evaluates this numerically:

-1/2 - I NIntegrate[DawsonF[Sqrt[u]]/(E^u - 1), {u, 0, Infinity}, 
    Method -> "DoubleExponential", WorkingPrecision -> 50]/Sqrt[Pi]
-1/2 - 0.93001253961059515359034795785857166233326206076173 I

Simultaneously computing a complete elliptic integral and its complement

2010-08-03T09:42:43Z

The complete elliptic integral of the first kind

$K(m)=\int_0^{\pi/2}\frac{\mathrm{d}t}{\sqrt{1-m\sin^2t}}$

is easily computed via the arithmetic-geometric mean iteration; to wit,

$K(m)=\frac{\pi}{2M(1,\sqrt{1-m})}$

where $M(a,b)$ is the arithmetic-geometric mean of $a$ and $b$ . With a little more trickery, the iteration can be hijacked to compute the complete elliptic integral of the second kind $E(m)$ as well.

In a number of applications, it happens that one needs both the values of $K(m)$ and its complement $K(1-m)$ (and sometimes similarly for $E(m)$ and $E(1-m)$ ).

My question is, apart from having to do an AGM iteration for each of $K(m)$ and $K(1-m)$ , is there an algorithm (maybe a modification of the basic AGM iteration) that simultaneously generates both $K(m)$ and its complement? I would also be interested in seeing also an extension of this algorithm, if one exists, for computing $E(m)$ as well (after which $E(1-m)$ is easily computed via Legendre's relation).

Answer by J. M. for The Riemann's Zeta Function represented as a continued fraction and a question of convergence.

2011-12-23T01:33:40Z

(Too long for a comment.)

There's a (somewhat) simpler (Eulerian) continued fraction:

$$\sum_{k=1}^\infty \frac1{k^s}=\sum_{k=1}^{\infty} \prod_{j=2}^k \left(1-\frac1{j}\right)^s=\cfrac1{1-\cfrac{\left(1-\frac12\right)^s}{1+\left(1-\frac12\right)^s-\cfrac{\left(1-\frac13\right)^s}{1+\left(1-\frac13\right)^s-\cfrac{\left(1-\frac14\right)^s}{1+\left(1-\frac14\right)^s-\cdots}}}}$$

but as you can see from comparing successive convergents of this continued fraction and the successive partial sums of the Dirichlet series, it's not terribly useful.

Also,

$$e^{-2z\,\mathrm{arcoth}(2k+1)}=\left(\frac{k}{k+1}\right)^z$$

so your CF could certainly be simplified a fair bit...

Answer by J. M. for Generalized trigonometric functions $Cos(n) v$ and $Sin(n) v$.

2011-12-17T02:23:46Z

(Too long for a comment.)

It's a bit older than your reference, but so-called "hypergoniometric functions" have been considered by Erik Lundberg in 1879. This article is a more recent discussion. Shelupsky and Burgoyne discuss similar generalizations. All ultimately consider this as the problem of inverting an appropriate generalization of the integral representations of arcsine and arccosine.

The $n=3$ case has been considered separately by A.C. Dixon; I had talked a bit about Dixon elliptic functions in this math.SE answer.

Answer by J. M. for Deriving Inverse of Hilbert Matrix

2010-11-28T13:00:01Z

It's a bit circuitous, but I'd like to point out this paper by Hitotumatu where he derives explicit expressions for the Cholesky triangle of a Hilbert matrix. From the expressions for the Cholesky triangle, you should be able to derive explicit expressions for the inverse (if $\mathbf A=\mathbf G\mathbf G^\top$ , then $\mathbf A^{-1}=\mathbf G^{-\top}\mathbf G^{-1}$ ).

Since the paper isn't that easily accessible, I'll include the main result here. If $\mathbf A$ is the Hilbert matrix, with the decomposition $\mathbf A=\mathbf L\mathbf D\mathbf L^\top$ with $\mathbf L$ unit lower triangular and $\mathbf D$ diagonal, then

$$\begin{align*}\ell_{j,k}&=\frac{(2k-1)\binom{2k-2}{k-1}\binom{2j-1}{j-k}}{(2j-1)\binom{2j-2}{j-1}}\\d_{k,k}&=\frac1{(2k-1)\binom{2k-2}{k-1}^2}\end{align*}$$

Answer by J. M. for interlacing roots/eigenvalues results and modern analogues

2011-11-22T03:08:07Z

A rather tenuous connection, but I'll throw it out and hope somebody can build a better answer from this (in short, too long for a comment):

It is well known (or it should be!) that one can use the Sturm sequence construction for a polynomial $p(x)$ with all roots real and its derivative to construct a symmetric tridiagonal companion matrix (that is, the symmetric tridiagonal matrix whose characteristic polynomial is $p(x)$). See Fiedler's paper for details.

One could then apply Cauchy's interlacing theorem on the tridiagonal matrix just constructed. As a matter of fact, the characteristic polynomials of successive leading submatrices of a tridiagonal form a Sturm sequence.

To expand on the comment I gave regarding orthogonal polynomials: consider the monic orthogonal polynomials $p_n(x)$ satisfying the difference equation

$$p_{n+1}(x)=(x-c_n)p_n(x)-d_np_{n-1}(x),\qquad p_{-1}(x)=0,p_0(x)=1$$

One can associate a symmetric tridiagonal matrix with this recursion, called a Jacobi matrix:

$$\begin{pmatrix}c_0&\sqrt{d_1}&&\\\sqrt{d_1}&\ddots&\ddots&\\&\ddots&&\sqrt{d_{n-1}}\\&&\sqrt{d_{n-1}}&c_{n-1}\end{pmatrix}$$

whose characteristic polynomial is $p_n(x)$. It can be seen that the characteristic polynomials of the leading $1\times1, 2\times 2, \dots$ submatrices are $p_1(x),p_2(x),\dots$ We see here the correspondence between the Sturm sequences for orthogonal polynomials and tridiagonal matrices.

Differential equation for a ratio of consecutive Bessel functions

2010-08-22T02:46:03Z

My attempts to search via Google seem to be failing, so I thought of asking here.

All the derivatives of the function

$r_n(z):=\frac{J_n(z)}{J_{n-1}(z)}$

where $J_n(z)$ is the Bessel function of the first kind are expressible in terms of $r_n(z)$ , for instance $\frac{\mathrm{d}}{\mathrm{d}z}r_n(z)=r_n(z)^2-\frac{2n-1}{z}r_n(z)+1$ . I've been trying to derive a (linear?) differential equation (hopefully just second-order) that might be satisfied by $r_n(z)$ , but my manipulative ability does not seem to be up to snuff.

Probably my problem can be resolved in two ways:

Are there any papers where generating functions/differential equations of ratios of Bessel functions have been studied? ; or
How can I derive a differential equation for $r_n(z)$ with the knowledge that all higher derivatives are expressible in terms of $r_n(z)$ ?

(On the other hand, the difference equation for $r_n(z)$ (and thus its continued fraction representation) is easily derived, so no problem for me there.)

I will be interested in any input. Thanks!

EDIT:

Per Pietro's request, I now tip my hand and reveal my reason for interest: I saw this paper many years ago on a neat method for computing the first few roots of the Bessel function of the first kind. Some time later, I came across J.F. Traub's "Iterative Methods for the Solution of Nonlinear Equations", where he shows the construction of iteration functions involving derivatives and can be constructed to have quadratic, cubic... convergence. (Newton's method is but the first member of this family). I also came across this short note by D.J. Hofsommer on how one might profitably exploit the methods derived by Traub if the function of interest satisfies a simple differential equation (Essentially, one just constructs the Newton correction $u=\frac{f(z)}{f^{\prime}(z)}$ , and the high-order iteration functions are merely a series in powers of $u$ ). That got me wondering on how one might recursively generate iteration functions with increasing order of convergence for the case of finding the roots of the Bessel function. (On another note, I was able to successfully use the ideas of Traub and Hofsommer for the generation of Gaussian quadrature rules, e.g. Legendre, Lobatto, Radau, and was hoping things might be just as successful for Bessel function root-finding).

nth-order generalizations of the arithmetic-geometric mean

2010-09-03T04:01:48Z

The arithmetic-geometric mean,

$a_{k+1}=\frac{a_k+b_k}{2} \quad b_{k+1}=\sqrt{a_k b_k}$

is one of the celebrated discoveries of Gauss, who found out that it is equivalent to computing a (complete) elliptic integral (which is a special case of the Gauss hypergeometric function ${}_2 F_1$).

I have been wondering if nth-order generalizations of the iteration,

$a_{k+1}=\frac{a_k+b_k}{n} \quad b_{k+1}=\sqrt[n]{a_k b_k}$

have ever been systematically studied. I've seen this paper by Borwein, but have had trouble searching for other papers. In particular, I'm interested if the coupled sequences also have a common limit, and if so, whether the limit is expressible as a hypergeometric function (or generalizations like those of Appell or Lauricella).

Another possible generalization I thought involves $n$ variables and makes use of the elementary symmetric polynomials. To use $n=4$ as an example:

$a_{k+1}=\frac{a_k+b_k+c_k+d_k}{4}$

$b_{k+1}=\sqrt{\frac{a_k b_k+a_k c_k+a_k d_k+b_k c_k+b_k d_k+c_k d_k}{3}}$

$c_{k+1}=\sqrt[3]{\frac{{a_k b_k c_k}+{a_k b_k d_k}+{a_k c_k d_k}+{b_k c_k d_k}}{2}}$

$d_{k+1}=\sqrt[4]{a_k b_k c_k d_k}$

Would these four sequences (and in general the $n$ sequences) tend to a common limit $F(a_0,b_0,c_0,d_0,\dots)$ like in the $n=2$ case, and if so, are they expressible in terms of known functions?

EDIT

Taking into account Darsh Ranjan's comments, I realized that what I should be looking at instead is the generalization whose denominators are binomial coefficients (thus, the general form $\sqrt[j]{\frac{e_j}{\binom{n}{j}}}$, for $j=1\dots n$ where $e_j$ is the jth elementary symmetric polynomial). The case $n=4$ now looks like

$a_{k+1}=\frac{a_k+b_k+c_k+d_k}{4}$

$b_{k+1}=\sqrt{\frac{a_k b_k+a_k c_k+a_k d_k+b_k c_k+b_k d_k+c_k d_k}{6}}$

$c_{k+1}=\sqrt[3]{\frac{{a_k b_k c_k}+{a_k b_k d_k}+{a_k c_k d_k}+{b_k c_k d_k}}{4}}$

$d_{k+1}=\sqrt[4]{a_k b_k c_k d_k}$

So, still the same question: is there a common limit, and if so, is the limit expressible in terms of known functions?

Answer by J. M. for Orthogonal polynomials/functions on the interval [0,1] but with same weight as Gegenbauer polynomials

2011-11-22T04:32:23Z

There is this paper and this paper which treat the special case of "half-range Chebyshev polynomials" (both kinds, corresponding to the weights $\dfrac1{\sqrt{1-x^2}}$ and $\sqrt{1-x^2}$ over $[0,1]$) to deal with Fourier expansions of nonperiodic functions. I have a feeling that half-range Gegenbauer polynomials have been treated before, and I'll try to see what I can dig up.

In the meantime, one can use the Stieltjes procedure to build up the recursion relations for these half range Gegenbauers. Letting

$$\langle f(x),g(x) \rangle^{(\alpha)}=\int_0^1 (1-t^2)^{\alpha-1/2} f(t)g(t)\mathrm dt$$

be the associated inner product, the Stieltjes procedure for generating monic orthogonal polynomials $\phi_k(x)$ uses the formulae

$$\begin{align*}b_k&=\frac{\langle x\phi_k(x),\phi_k(x)\rangle^{(\alpha)}}{\langle\phi_k(x),\phi_k(x)\rangle^{(\alpha)}}\\ c_k&=\frac{\langle\phi_k(x),\phi_k(x)\rangle^{(\alpha)}}{\langle\phi_{k-1}(x),\phi_{k-1}(x)\rangle^{(\alpha)}}\end{align*}$$

to give the coefficients $b_k,c_k$ for the recursion relation

$$\phi_{k+1}(x)=(x-b_k)\phi_k(x)-c_k\phi_{k-1}(x)$$

Here, the result

$$\int_0^1 (1-t^2)^{\alpha-1/2}t^k \mathrm dt=\frac{\Gamma\left(\frac{1+k}{2}\right)\Gamma\left(\alpha+\frac12\right)}{2\Gamma\left(\alpha+\frac{k}{2}+1\right)}$$

is useful.

I might as well throw this in. There is an algorithm due to Chebyshev (1859) for determining recursion coefficients from the moments. I've already talked about the algorithm here, so I shall not repeat myself. Instead, I'll reproduce the Mathematica routine I gave in that answer:

chebAlgo[mom_?VectorQ, prec_: MachinePrecision] := 
 Module[{n = Quotient[Length[mom], 2], si = mom, ak, bk, np, sp, s, v},
  np = Precision[mom]; If[np === Infinity, np = prec];
  ak[1] = mom[[2]]/First[mom]; bk[1] = First[mom];
  sp = PadRight[{First[mom]}, 2 n - 1];
  Do[
   sp[[k - 1]] = si[[k - 1]];
   Do[
    v = sp[[j]];
    sp[[j]] = s = si[[j]];
    si[[j]] = si[[j + 1]] - ak[k - 1] s - bk[k - 1] v;
    , {j, k, 2 n - k + 1}];
   ak[k] = si[[k + 1]]/si[[k]] - sp[[k]]/sp[[k - 1]];
   bk[k] = si[[k]]/sp[[k - 1]];
   , {k, 2, n}];
  N[{Table[ak[k], {k, n}], Table[bk[k], {k, n}]}, np]
  ]

Here for instance is how to use chebAlgo[] to generate recursion coefficients for the monic half-range Chebyshev polynomials of the first kind:

With[{a = 0}, chebAlgo[Table[Gamma[(k + 1)/2] Gamma[a + 1/2]/Gamma[a + k/2 + 1],
                             {k, 0, 10}]/2, Infinity]] // FullSimplify

Answer by J. M. for Sum of products of exponentials and polynomials

2011-11-22T06:24:19Z

You really can't do better than the polylogarithm and Lerch's transcendent as closed forms.

For what it's worth, there are the explicit representations

$$\begin{align*}\mathrm{Li}_{-n}(z)&=\frac1{(1-z)^{n+1}} \sum_{m=1}^n \left(\sum_{k=1}^m (-1)^{k+1} \binom{n+1}{k-1}(m-k+1)^n\right)z^m\\\Phi(z,-n,a)&=a^n+\sum_{j=0}^n \binom{n}{j} \mathrm{Li}_{-j}(z) a^{n-j}\end{align*}$$

all valid for $n$ a positive integer. Daunting, no?

If you can't accept these closed forms, then you're better off staying with the sum representation that you have...

On a polynomial related to the Legendre function of the second kind

2011-09-04T02:27:04Z

The Legendre function of the second kind, $Q_n(z)$, along with the usual Legendre polynomial $P_n(z)$, are the two linearly independent solutions of the Legendre differential equation.

$Q_n(z)$ can be expressed in the following form:

$$Q_n(z)=P_n(z)\mathrm{artanh}\,z-W_{n-1}(z)$$

where $W_{n-1}(z)$ can be expressed either as

$$W_{n-1}(z)=\sum_{k=1}^n \frac{P_{k-1}(z) P_{n-k}(z)}{k}$$

or as

$$W_{n-1}(z)=\sum_{k=0}^{n-1} \frac{(H_n-H_k)(n+k)!}{2^k (n-k)! (k!)^2} (z-1)^k$$

where $H_k$ is the $k$-th harmonic number, $H_k=\sum\limits_{j=1}^k \frac1{j}$.

My questions:

Mathematica returns a rather complicated expression for $W_{n-1}(z)$ involving the unknown solution of a certain recurrence (i.e. DifferenceRoot[]). Is there possibly a simpler form for this polynomial?
Might there be a (hopefully simple) $n$-term recurrence that generates these polynomials?

Addendum:

After staring long and hard at Pietro's answer, I feel now that my second question was sorta kinda dumb; I already knew that both Legendre functions satisfied the same difference equation, so it stands to reason that a linear combination of them should also be a solution to that recurrence.

I now would like to expand my first question a bit: is it possible to express $W_n(z)$ as a single hypergeometric function (e.g. ${}_p F_q$ or some of the fancy multivariate ones), perhaps with one of the parameters being a negative integer? For instance, $P_n(z)$ is expressible as a Gaussian hypergeometric function ${}_2 F_1$ with one of the numerator parameters being a negative integer. Might there be something similar for the $W_n$?

I would also like to consider an additional question: are the $W_n$ orthogonal polynomials with respect to some weight function $\omega(x)$ and an associated support interval $(a,b)$? That is, if

$$\int_a^b\omega(t)W_j(t)W_k(t)\mathrm dt=0,\qquad j\neq k$$

for some $\omega(x)$ and some interval $(a,b)$, what is this weight function and its support interval?

Answer by J. M. for roots of polynomial with matrix coefficients

2011-08-21T13:08:58Z

Yes, there a number of numerical methods for finding the solvents of a polynomial with matrix coefficients. Dennis, Traub, and Weber in this article give some of the relevant theory, as well as some algorithms for finding so-called "dominant solvents" for your matrix polynomial (see also their earlier article). This article presents the use of Newton's method for solving matrix polynomials.

As a tiny aside on the "companion" formulation mentioned by Federico: make sure that whatever software you're using is the version that uses the QZ algorithm, without the preliminary application of the inverse of the leading coefficient to "monicize" the matrix polynomial. It can happen that the leading coefficient is an ill-conditioned matrix, and the "monicization" leads to a degradation of accuracy in the solvents returned.

How many curves can fit on a sphere without intersecting?

2011-08-17T03:07:14Z

(This is more or less a generalized and more abstract version of this m.SE question.)

Let's recall two facts first:

Any (arc-length parametrized) space curve is uniquely determined (up to rigid motions) by its curvature $\kappa(s)$ and its torsion $\tau(s)$.

The curvature and torsion of a unit-spherical curve satisfies the equation $$\left(\tau\frac{\mathrm d}{\mathrm ds}\frac1{\kappa}\right)^2+\frac1{\kappa^2}=1$$

In my m.SE question, the unique properties of the loxodrome and the existence of the Mercator projection allowed the determination of when $n$ loxodromes can be placed on a sphere without intersecting. I am now considering the following question:

What conditions should the curvature and torsion of a spherical curve satisfy so that $n$ copies of it can be placed on a sphere without intersecting?

Alternatively, any spherical curve can also be defined as parametric equations of the longitude $\theta=\theta(t)$ and co-latitude $\varphi=\varphi(t)$, or in Cartesian form (for a unit sphere):

$$\begin{align*}x&=\cos\,\theta(t)\sin\,\varphi(t)\\y&=\sin\,\theta(t)\sin\,\varphi(t)\\z&=\cos\,\varphi(t)\end{align*}$$

so we can also ask

What conditions should $\theta(t)$ and $\varphi(t)$ satisfy so that $n$ copies of it can be placed on a sphere without intersecting?

If there are better ways to characterize the conditions than these, I'd be interested in hearing about them.

For convenience, let's consider only the unit sphere.

One theorem(?) I have along these lines is

If a spherical curve fits within a hemisphere without touching a great circle, then two copies of the curve can be fit on the sphere without intersecting.

but I haven't figured out how to generalize this.

Answer by J. M. for do numerical integration with fixed abscissas

2011-08-05T05:34:46Z

I have settled this question here; in brief, a FORTRAN implementation of algorithms for solving this problem have been published in the Collected Algorithms of the ACM.

At the risk of being redundant, I'd like to mention here some other things I mentioned in the (now expanded) answer of mine at that other site for completeness' sake.

Robert's warning of the Runge phenomenon happening is a good one, and it does happen if your abscissas are perversely distributed (relatedly, the underlying Vandermonde matrix is ill-conditioned); the equispaced case being among the worst-behaved point distributions. Abscissas that are "nicely distributed" (e.g. Legendre, Chebyshev, or any other point distribution which "clusters" near the ends) will generally ensure that you have a quadrature rule that behaves tamely even for large numbers of points. (As an aside, the chebfun project hinges on functions being nicely approximated by interpolating polynomials with abscissas that are transformed Chebyshev polynomial roots/extrema.)

Lastly, whatever you finally settle with, you will want to perform the sanity check of ensuring that all the weights of your quadrature rule are of the same sign; any change of sign in the weights can lead to subtractive cancellation when you employ the quadrature rule, and you wouldn't want that...

Answer by J. M. for Valid to use all Wynn-extrapolated values?

2011-08-05T05:54:59Z

There really isn't a theorem you can use, since most of the applicable theory is for sequences with known asymptotic behavior. For "in the wild" sequences, you can do no better than to check that the results of the $\varepsilon$-algorithm remain sensible as you proceed. As a general rule however, for the recursion

$$\varepsilon_{k+1}^{(n)}=\varepsilon_{k-1}^{(n+1)}+\frac1{\varepsilon_{k}^{(n+1)}-\varepsilon_{k}^{(n)}}$$

the so-called "diagonal approximations" (borrowing from the theory of Padé approximants) $\varepsilon_m^{(0)}$ or $\varepsilon_m^{(1)}$, depending on which diagonal of the table you're already in, would be the most accurate approximations for the limit of your sequence. Still, for these algorithms, the human eye is a better judge of convergence than the computer, so have your Wynn routine print out the entire array if it can.

For more practical advice, have a look at E.J. Weniger's survey paper.

Answer by J. M. for Numerical Beta Function

2011-08-05T05:38:33Z

There's nothing more straightforward than using the gamma function relationship, I believe (perhaps using Lanczos's approximation to compute the gamma functions). Of course, for integer arguments, the product representation is much faster. You could probably also consider the special case when one of the beta function's arguments is a semi-integer, and accordingly construct the appropriate product representation (involving $\sqrt{\pi}$).

Finding a recursion for a sum of Legendre polynomials

2010-08-20T17:45:42Z

The polynomial

$a_n(x):=P_n(x)-\frac{n-1}{n}P_{n-2}(x)$

where $P_n(x)$ is a Legendre polynomial came up while I was investigating methods for estimating the error in Gaussian quadrature.

I am wondering if a recursion relation (homogeneous/inhomogeneous, three terms or more) might be found for the $a_n(x)$. Of course I can obviously just generate them directly after performing the recursion for Legendre polynomials, but I'm wondering if there might be a way to generate this sequence of polynomials from the first few members, for instance $a_1(x)=x$ and $a_{2}(x)=\frac{3}{2}x^2-1$...

I have heard of Sister Celine's method; my understanding of it however is that one must first have a representation for your polynomials of interest in terms of a hypergeometric function. So an alternative question might be: can $a_n(x)$ be represented as a single hypergeometric function, most probably with one of the numerator parameters being -n?

Estimating the spectral radius of a matrix, noniteratively

2010-08-13T06:31:10Z

Morris Marden's "Geometry of Polynomials" displays a number of formulae that allow one to estimate bounds on the largest root of a polynomial that do not involve actual rootfinding. Having been inspired by this, and since this particular problem crops up in one of the things I'm working on, I was wondering if one could get good estimates of the spectral radius of a general dense n-by-n matrix A that has been previously processed as follows:

a similarity transformation to upper Hessenberg form ( $A=QHQ^T$ , $Q$ orthogonal and $H$ upper Hessenberg); and
subtracting the identity multiplied by the mean of the eigenvalues from $H$ ( $H'=H-\frac{trace(H)}{n}I$ , this corresponds geometrically to centering the eigenvalues around the origin of the complex plane).

As much as possible, I am trying to avoid having to resort to an eigenvalue method (e.g. QR (too much effort!), power method (the power method can misbehave when there is more than one eigenvalue whose modulus is equal to the spectral radius)) since I only need a quick 2-3 digit approximation of the spectral radius. I have considered actually expanding $H'$ to its characteristic polynomial (equivalently, a similarity transformation of $H'$ to a Frobenius companion matrix) so that the formulae listed in Marden can apply, but after reading Wilkinson's wonderful book "The Algebraic Eigenvalue Problem" where he details how unstable the computation of coefficients of the characteristic polynomial can get from a matrix, I suppose that idea is shot.

My other naïve attempts at estimating the spectral radius include using $||H'||_\infty$ as an estimate, and deriving rough bounds using Gerschgorin's theorem; the problem I've seen is that both attempts tend to overestimate the spectral radius by a significant factor.

Is there a way to estimate the spectral radius more cheaply and noniteratively than actually computing eigenvalues?

Answer by J. M. for Efficiently computing a matrix's induced p-norm

2010-09-18T01:06:18Z

Nicholas Higham gives an algorithm for estimating the Hölder $p$-norm of a matrix with the estimate being within a factor of $n^{1-1/p} \|\mathbf{A}\|_p$ ; maybe you can somehow adapt this approach to your needs?

(added 5/13/2011)

I posted a Mathematica translation of Higham's original MATLAB code here.

What's the difference between a Riemann theta and a Siegel theta function?

2011-05-08T06:55:13Z

One of the things I'm working on has required me to look into the literature of multidimensional theta functions, and I've gotten a bit confused on a few naming details.

A look at the DLMF says that "the" multidimensional theta function is the Riemann theta function,

$$\mathop{\Theta}\left(\mathbf{z}\mid\boldsymbol{{\Omega}}\right)=\sum_{{\mathbf{n}\in{\mathbb Z}^g}}e^{{2\pi i\left(\frac{1}{2}\mathbf{n}\cdot\boldsymbol{{\Omega}}\cdot\mathbf{n}+\mathbf{n}\cdot\mathbf{z}\right)}}$$

with $\mathbf{z}$ a complex vector and $\mathbf \Omega$ a complex symmetric matrix with symmetric positive definite imaginary part. I find that there is a nice implementation for numerically computing this function, but it is only available in Maple and Java. Since I work in Mathematica, I decided to see if there was already something in Mathematica before going through the trouble of translating code.

Peering at the docs for Mathematica netted the Siegel theta function. However, looking at how the Siegel theta function was defined, this and the Riemann theta function seem to be the same thing!

Not having Maple to check if the results from their implementation of Riemann theta and the results from Mathematica's implementation of Siegel theta agree, and not being able to access the (older) references to these functions, I now wish to ask: is there no difference whatsoever between these two except in name (and if there are, how are these two different)? Did Riemann and Siegel independently study the exact same multidimensional theta function? How come it seems that there is no reference to these two being synonymous?

Answer by J. M. for Numerical Differentiation. What is the best method?

2011-05-08T15:33:31Z

If your function is badly behaved (e.g. noisy, very oscillatory), no method will perform properly (differentiation is numerically very unstable). That being said, for "nice functions", I have good experience with polynomial (Richardson) extrapolation methods. This paper and this paper give hints on how you might write your own implementation. I will note that this is the method implemented in the NAG numerical libraries (with of course a few wrinkles of their own).

There are two possible alternatives if for some reason you don't want to use Richardsonian methods. One is to use Cauchy's differentiation formula:

$$f^\prime(x)=\frac1{2\pi i}\oint_\gamma \frac{f(t)}{(t-x)^2}\mathrm dt$$

where it is up to you to choose a suitable counterclockwise contour $\gamma$ (a circle is customary); the other is to use the "Lanczos derivative":

$$f^\prime(x)=\lim_{h\to 0}\frac{3}{2h^3}\int_{-h}^h t\;f(x+t)\mathrm dt$$

where you either will have to experiment with an appropriate step size $h$, or use some extrapolative procedure.

You will have to experiment with your computing environment to choose.

Some Dirichlet series questions.

2010-12-28T06:31:24Z

I asked this question on m.SE in an attempt to find out the right words to say for these questions I am about to ask.

In his great answer, Matthew Emerton explained that (cuspidal) automorphic L-functions correspond to the Dirichlet series with "nice" properties like having a reflection equation, a meromorphic continuation to the entire complex plane and a suitable analogue of the Riemann hypothesis. This leads me to my first question:

1) Are there Dirichlet series that cannot be classified as (cuspidal) automorphic L-functions, yet still possess a critical line of nontrivial zeroes?

Now, I come to the question I had been meaning to ask here. It is known that functions like Riemann $\zeta$ and the Ramanujan Dirichlet series admit a "Riemann-Siegel" decomposition; that is, letting $\sigma$ denote the position of the "critical line" of the Dirichlet series $g(s)$, they can be expressed as

$$g(\sigma+it)=z(t)\exp(-i\vartheta(t))$$

where $z(t)$ and $\vartheta(t)$ are "Riemann-Siegel" functions corresponding to the Dirichlet series $g(s)$. The value of $z(t)$ is that it eases the task of finding nontrivial zeroes of the corresponding Dirichlet series (essentially helping to verify its corresponding "hypothesis"). My question now is

2) Do all (cuspidal) automorphic L-functions have a "Riemann-Siegel" decomposition? If not, what restrictions are there for them to possess such a decomposition?

My motivation is more of curiosity than anything else. Hopefully this is not too elementary a question!

Answer by J. M. for What would you want to see at the Museum of Mathematics?

2010-12-25T16:28:51Z

Sculptures of surfaces would be lovely.

Answer by J. M. for Second order Taylor expansion to solve system of equations

2010-12-22T15:56:43Z

I haven't gotten around to downloading and reading it (and I'm wondering how I missed this when I was searching for results related to Halley's method), but apparently a multivariate version of the Halley iteration has already been developed decades ago. Maybe this might be of use.

Answer by J. M. for Why not evaluate integrals using ODE-solvers?

2010-11-30T09:42:14Z

Here's my take on the matter: the difference of philosophy between quadrature routines and ODE solving routines, I believe, is this:

Extrapolation is riskier than interpolation.

Remember that numerical quadrature routines all boil down to approximating your presumably more complicated integrand within the interval of integration with something easier to integrate exactly, and then integrating that. For instance, Newton-Cotes (and in essence, Romberg as well) constructs an interpolating polynomial from your integrand with equispaced abscissas, and integrating that. For Gaussian or Clenshaw-Curtis quadrature, it is equivalent to interpolating your function at "specially spaced" abscissas (Legendre polynomial roots in the former, and Chebyshev polynomial roots in the latter) that have better convergence in the limit. In effect, we run under the assumption that the interpolating function behaves very similarly to the actual integrand within the interval of interest that a sufficient amount of sampling within the integration interval should be enough to capture the behavior of your integrand, and thus give a result hopefully close to the actual value of the integral.

In contrast, remember that ODE solvers usually only have initial values to start with. The reason for building in a lot of machinery in current ODE solvers, whether Runge-Kutta, Bulirsch-Stoer, Adams/Gear multistep, or some of the fancier modern techniques, is that extrapolation is inherently unstable. Knowing how the solution looks like in the beginning gives no guarantee how it will behave as the ODE solver marches on; the solution may well be violently oscillatory, or decaying quite fast (so-called "stiff" problems). Thus, there is quite a fair amount of code inscribed in modern ODE solvers for checking how reasonable are the step-sizes being taken, and other such fail-safes.

As I did mention in some previous comments, some ODE solving methods are equivalent to quadrature methods when applied to the initial-value problem $y^{\prime}=f(x)$: using classical Runge-Kutta for quadrature is equivalent to performing Simpson's rule, for instance.

The point is that ODE solvers tend to be more careful ("tiptoeing", if you will) and thus more effort-intensive than numerical quadrature routines because they make no assumptions on how your integrand behaves. On that note, I will say that you should know that there are integrands (and corresponding intervals) where using an ODE solver might make more sense than using a numerical quadrature routine. One instance that comes to mind: if you know (through graphing, for instance) that your integrand has crazy behavior in a relatively tiny interval within the interval of integration, while the sampling done by a numerical quadrature routine might miss such features (or take a long time to notice them), an ODE solver will be careful to begin with and is less likely to miss the crazy behavior, and will shrink step-sizes as appropriate until it has gone past that interval.

Answer by J. M. for Fast trace of inverse of a square matrix

2010-11-19T12:35:07Z

Given that the poster has specified that his matrix is symmetric, I offer a general solution and a special case:

Eigendecomposition actually becomes more attractive here: the bulk of the work is in reducing the symmetric matrix to tridiagonal form, and finding the eigenvalues of a tridiagonal matrix is an O(n) process. Assuming that the symmetric matrix is nonsingular, summing the reciprocals of the eigenvalues nets you the trace of the inverse.
If the matrix is positive definite as well, first perform a Cholesky decomposition. Then there are methods for generating the diagonal elements of the inverse.

Answer by J. M. for Most memorable titles

2010-10-31T14:46:55Z

The AKS paper Primes is in P is a pretty memorable title for me.

Answer by J. M. for Square root of non-positive definite matrix

2010-10-28T09:08:21Z

Personally, I'm still a bit torn while giving this prescription. I am of the opinion that for truly reliable computation, one should use the singular value decomposition whenever one can tolerate the performance penalty; it truly is a very stable algorithm that can also be used to compute a bunch of important diagnostic quantities. (As an aside, if it is truly the square root you want, you could exploit the fact that having a singular value decomposition of $\mathbf{A}$ is effectively equivalent to having an eigendecomposition of either of $\mathbf{A}^T \mathbf{A}$ or $\mathbf{A}\mathbf{A}^T$)

Having said this, if you're still dead-set on using Cholesky on a positive semidefinite matrix, while in exact arithmetic you are supposed to encounter a zero, what might actually happen with a (large enough) matrix with inexact entries is that your Cholesky routine encounters a tiny quantity not detected as zero, and the routine happily carries over this tiny quantity to divide other entries with. Disaster!

The cure is that one does a symmetric pivoting $\mathbf{A}\to\mathbf{P}\mathbf{A}\mathbf{P}^T$ (where $\mathbf{P}$ is a permutation matrix), which reorders the diagonal entries of $\mathbf{A}$ (no off-diagonal entries are moved into the diagonal). This has the effect that the (near-)zero quantities are not encountered until one already has proceeded through the $r$-th iteration of the main loop in the Cholesky decomposition, where $r$ is the (perceived) rank of the matrix (this depends on what tolerance you have set, or the value such that anything whose absolute value is lower than it is effectively treated as zero).

I have only given a brief sketch, since what I really should be doing is to point out this paper by Nick Higham, which discusses the nuances of Cholesky decomposition with symmetric pivoting for symmetric positive semidefinite matrices.

Answer by J. M. for How do eigenvectors and eigenvalues change when we remove a row/column pair of a matrix?

2010-10-25T13:36:05Z

The word you're looking for is downdating, and I cannot do better than to point out these two survey papers, and this article. I will also have to make the reminder that it makes better numerical sense to compute the singular values of $\mathbf{D}$ rather than the eigenvalues of $\mathbf{D}^T \mathbf{D}$.

Comment by J. M.

2013-05-30T19:04:25Z

(of course, since all I put in my answer there were links to related papers, it can't really be said to be "better"...)

Comment by J. M.

2013-05-30T19:03:06Z

As you have linked to the paper by Taussky and Todd, you might want to mention www-math.mit.edu/~edelman/homepage/papers/kac.pdf as well...

Comment by J. M.

2011-12-23T01:50:10Z

That solution gives a complex root, not a real one.

Comment by J. M.

2011-12-23T01:36:56Z

I commented in his math.SE question that he try van Vleck, but I guess he wanted somebody to do it for him... :) Maybe the parabola theorem might yield something more useful, but I'm far away from my refs.

Comment by J. M.

2011-12-23T01:17:54Z

@Noam: it's somewhat standard in CF literature. I'm told it's originally Gauss's.

Comment by J. M.

2011-12-18T10:47:19Z

The article @Doug linked to is also at oberlin.edu/faculty/jcalcut/gausspi.pdf ; see also oberlin.edu/faculty/jcalcut/arctan.pdf .

Comment by J. M.

2011-12-16T14:32:34Z

uni-graz.at/~gronau/monthly230-237.pdf is a newer paper discussing the analytic continuation of the Theodorus spiral.

Comment by J. M.

2011-12-16T12:31:42Z

Mathematica and Maple sidestep the problem of "symbolic expression" of roots by using a data structure containing the minimal integer-coefficient polynomial that corresponds to those roots, and an index. For actual general expressions for roots of polynomials with degree $\geq 5$ , you need theta functions. Which is a pretty deep rabbit hole...

Comment by J. M.

2011-12-15T01:51:46Z

Eep, sorry for not checking that first DOI! Here's the first one: springerlink.com/content/u528652t80v96515

Comment by J. M.

2011-12-14T13:28:33Z

Hey, I don't know if you've already found an answer to this, but it looks to me that the results of dx.doi.org/10.1007/BFb0120751 and dx.doi.org/10.1137/0131057 might be applicable to your problem.

Comment by J. M.

2011-12-10T11:33:17Z

For cycle detection, either of Floyd's or Brent's algorithms should be suitable.

Comment by J. M.

2011-12-10T04:26:10Z

ren, you might want to also try asking this question on math.stackexchange.com ; one of the active answerers there is experienced with using Mellin transforms for integrals of this sort. Just make sure to mention that you've already asked here, and link to this question.

Comment by J. M.

2011-12-09T04:22:36Z

It does, but only after you build-in symmetric pivoting. Even then, it's easy to choke $LDL^\top$; try $\begin{pmatrix}0&1\\1&0\end{pmatrix}$.

Comment by J. M.

2011-12-09T03:57:14Z

@Will: in the case of the elliptic integral of the third kind $\Pi(n;\phi\mid m)$, depending on the characteristic $n$, the computations can either be rather simple or a bit complicated. See dx.doi.org/10.1007/BF02165405 and dx.doi.org/10.1007/BF02167558 for instance. There is now the algorithm due to Carlson: dx.doi.org/10.1007/BF02198293 , which is slightly less complicated, but I haven't done a comparison of this versus plain quadrature with, say, Gauss-Kronrod.

Comment by J. M.

2011-12-09T03:43:40Z

@aukm: Yes, I do believe you can economize that series you already have; alternatively, there is also the option of constructing the Padé approximant directly from your series. But having done so, I still think polishing the results of either of those with Newton-Raphson would be the best you can make of this.