User j. m. - MathOverflow most recent 30 from http://mathoverflow.net 2013-06-20T07:30:43Z http://mathoverflow.net/feeds/user/7934 http://www.creativecommons.org/licenses/by-nc/2.5/rdf http://mathoverflow.net/questions/131810/is-it-possible-to-generalize-functions-like-xy-ln-x-sin-x-arctan-x-to-su/132357#132357 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? J. M. 2013-05-30T18:53:44Z 2013-05-30T18:53:44Z <p>(Too long for a comment.)</p> <p>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</p> <p><code>$$f(z)=\frac1{2\pi i}\oint_\gamma f(t) (t-z)^{-1}\mathrm dt$$</code></p> <p>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.</p> http://mathoverflow.net/questions/104082/evaluating-a-limit-similar-to-the-euler-constant/104168#104168 Answer by J. M. for Evaluating a limit similar to the Euler constant J. M. 2012-08-07T01:19:05Z 2012-08-07T01:19:05Z <p>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</p> <p><code>$$-\frac{i}{2}\sum_{k=1}^\infty \frac1{ik+k^{3/2}}$$</code></p> <p>Using Laplace transform techniques, your sum can be transformed into the integral</p> <p><code>$$-\frac12-\frac{i}{\sqrt\pi}\int_0^\infty \frac{F(\sqrt{u})}{\exp\,u-1}\mathrm du$$</code></p> <p>where $F(z)$ is <a href="http://dlmf.nist.gov/7.2.E5" rel="nofollow">Dawson's integral</a>.</p> <p>I don't know of a closed form for this integral, but <em>Mathematica</em> easily evaluates this numerically:</p> <pre><code>-1/2 - I NIntegrate[DawsonF[Sqrt[u]]/(E^u - 1), {u, 0, Infinity}, Method -&gt; "DoubleExponential", WorkingPrecision -&gt; 50]/Sqrt[Pi] -1/2 - 0.93001253961059515359034795785857166233326206076173 I </code></pre> http://mathoverflow.net/questions/34363/simultaneously-computing-a-complete-elliptic-integral-and-its-complement Simultaneously computing a complete elliptic integral and its complement J. M. 2010-08-03T09:42:43Z 2012-04-30T11:17:39Z <p>The complete elliptic integral of the first kind</p> <p><code>$K(m)=\int_0^{\pi/2}\frac{\mathrm{d}t}{\sqrt{1-m\sin^2t}}$</code></p> <p>is easily computed via the arithmetic-geometric mean iteration; to wit,</p> <p><code>$K(m)=\frac{\pi}{2M(1,\sqrt{1-m})}$</code></p> <p>where <code>$M(a,b)$</code> is the arithmetic-geometric mean of <code>$a$</code> and <code>$b$</code>. With a little more trickery, the iteration can be hijacked to compute the complete elliptic integral of the second kind <code>$E(m)$</code> as well.</p> <p>In a number of applications, it happens that one needs both the values of <code>$K(m)$</code> and its complement <code>$K(1-m)$</code> (and sometimes similarly for <code>$E(m)$</code> and <code>$E(1-m)$</code>).</p> <p>My question is, apart from having to do an AGM iteration for each of <code>$K(m)$</code> and <code>$K(1-m)$</code>, is there an algorithm (maybe a modification of the basic AGM iteration) that simultaneously generates both <code>$K(m)$</code> and its complement? I would also be interested in seeing also an extension of this algorithm, if one exists, for computing <code>$E(m)$</code> as well (after which <code>$E(1-m)$</code> is easily computed via Legendre's relation).</p> http://mathoverflow.net/questions/84108/the-riemanns-zeta-function-represented-as-a-continued-fraction-and-a-question-of/84132#84132 Answer by J. M. for The Riemann's Zeta Function represented as a continued fraction and a question of convergence. J. M. 2011-12-23T01:33:40Z 2011-12-23T11:44:10Z <p>(Too long for a comment.)</p> <p>There's a (somewhat) simpler (Eulerian) continued fraction:</p> <p><code>$$\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}}}}$$</code></p> <p>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.</p> <p>Also,</p> <p><code>$$e^{-2z\,\mathrm{arcoth}(2k+1)}=\left(\frac{k}{k+1}\right)^z$$</code></p> <p>so your CF could certainly be simplified a fair bit...</p> http://mathoverflow.net/questions/83677/generalized-trigonometric-functions-cosn-v-and-sinn-v/83684#83684 Answer by J. M. for Generalized trigonometric functions $Cos(n) v$ and $Sin(n) v$. J. M. 2011-12-17T02:23:46Z 2011-12-17T02:23:46Z <p>(Too long for a comment.)</p> <p>It's a bit older than your reference, but so-called "hypergoniometric functions" have been considered by <a href="http://www.maths.lth.se/matematiklu/personal/jaak/hypergf.ps" rel="nofollow">Erik Lundberg in 1879</a>. <a href="http://www.jstor.org/pss/2695794" rel="nofollow">This article</a> is a more recent discussion. <a href="http://www.jstor.org/pss/2309789" rel="nofollow">Shelupsky</a> and <a href="http://dx.doi.org/10.1090/S0025-5718-1964-0164066-5" rel="nofollow">Burgoyne</a> discuss similar generalizations. All ultimately consider this as the problem of inverting an appropriate generalization of the integral representations of arcsine and arccosine.</p> <p>The $n=3$ case has been considered separately by <a href="http://gdz.sub.uni-goettingen.de/de/dms/load/img/?PPN=PPN600494829_0024&amp;PHYSID=PHYS_0179" rel="nofollow">A.C. Dixon</a>; I had talked a bit about Dixon elliptic functions in <a href="http://math.stackexchange.com/a/36506" rel="nofollow">this math.SE answer</a>.</p> http://mathoverflow.net/questions/47561/deriving-inverse-of-hilbert-matrix/47583#47583 Answer by J. M. for Deriving Inverse of Hilbert Matrix J. M. 2010-11-28T13:00:01Z 2011-12-07T01:52:37Z <p>It's a bit circuitous, but I'd like to point out this <a href="http://dx.doi.org/10.1007/BF03167904" rel="nofollow">paper</a> 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 <code>$\mathbf A=\mathbf G\mathbf G^\top$</code>, then <code>$\mathbf A^{-1}=\mathbf G^{-\top}\mathbf G^{-1}$</code>).</p> <hr> <p>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</p> <p><code>\begin{align*}\ell_{j,k}&amp;=\frac{(2k-1)\binom{2k-2}{k-1}\binom{2j-1}{j-k}}{(2j-1)\binom{2j-2}{j-1}}\\d_{k,k}&amp;=\frac1{(2k-1)\binom{2k-2}{k-1}^2}\end{align*}</code></p> http://mathoverflow.net/questions/81388/interlacing-roots-eigenvalues-results-and-modern-analogues/81571#81571 Answer by J. M. for interlacing roots/eigenvalues results and modern analogues J. M. 2011-11-22T03:08:07Z 2011-12-06T16:38:39Z <p>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):</p> <p>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 <a href="http://dx.doi.org/10.1016/0024-3795%2890%2990323-5" rel="nofollow">Fiedler's paper</a> for details.</p> <p>One could then apply <a href="http://books.google.com/books?id=mWin91cLndsC&amp;pg=PA202" rel="nofollow">Cauchy's interlacing theorem</a> on the tridiagonal matrix just constructed. As a matter of fact, <a href="http://books.google.com/books?id=mWin91cLndsC&amp;pg=PA141" rel="nofollow">the characteristic polynomials of successive leading submatrices of a tridiagonal form a Sturm sequence</a>.</p> <hr> <p>To expand on the comment I gave regarding orthogonal polynomials: consider the <em>monic</em> orthogonal polynomials $p_n(x)$ satisfying the difference equation</p> <p>$$p_{n+1}(x)=(x-c_n)p_n(x)-d_np_{n-1}(x),\qquad p_{-1}(x)=0,p_0(x)=1$$</p> <p>One can associate a symmetric tridiagonal matrix with this recursion, called a Jacobi matrix:</p> <p><code>$$\begin{pmatrix}c_0&amp;\sqrt{d_1}&amp;&amp;\\\sqrt{d_1}&amp;\ddots&amp;\ddots&amp;\\&amp;\ddots&amp;&amp;\sqrt{d_{n-1}}\\&amp;&amp;\sqrt{d_{n-1}}&amp;c_{n-1}\end{pmatrix}$$</code></p> <p>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.</p> http://mathoverflow.net/questions/36350/differential-equation-for-a-ratio-of-consecutive-bessel-functions Differential equation for a ratio of consecutive Bessel functions J. M. 2010-08-22T02:46:03Z 2011-12-02T02:24:51Z <p>My attempts to search via Google seem to be failing, so I thought of asking here.</p> <p>All the derivatives of the function</p> <p><code>$r_n(z):=\frac{J_n(z)}{J_{n-1}(z)}$</code></p> <p>where <code>$J_n(z)$</code> is the Bessel function of the first kind are expressible in terms of <code>$r_n(z)$</code>, for instance <code>$\frac{\mathrm{d}}{\mathrm{d}z}r_n(z)=r_n(z)^2-\frac{2n-1}{z}r_n(z)+1$</code> . I've been trying to derive a (linear?) differential equation (hopefully just second-order) that might be satisfied by <code>$r_n(z)$</code>, but my manipulative ability does not seem to be up to snuff.</p> <p>Probably my problem can be resolved in two ways:</p> <ol> <li>Are there any papers where generating functions/differential equations of ratios of Bessel functions have been studied? ; or</li> <li>How can I derive a differential equation for <code>$r_n(z)$</code> with the knowledge that all higher derivatives are expressible in terms of <code>$r_n(z)$</code>?</li> </ol> <p>(On the other hand, the <em>difference</em> equation for <code>$r_n(z)$</code> (and thus its continued fraction representation) is easily derived, so no problem for me there.)</p> <p>I will be interested in any input. Thanks!</p> <p><strong>EDIT:</strong></p> <p>Per Pietro's request, I now tip my hand and reveal my reason for interest: I saw <a href="http://dx.doi.org/10.1016/S0010-4655%2898%2900193-3" rel="nofollow">this paper</a> 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 <a href="http://books.google.com/books?id=se3YdgFgz4YC&amp;pg=PA78" rel="nofollow">construction of iteration functions involving derivatives</a> 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 <a href="http://www.ams.org/journals/mcom/1958-12-061/S0025-5718-1958-0099752-6/S0025-5718-1958-0099752-6.pdf" rel="nofollow">short note</a> 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 <code>$u=\frac{f(z)}{f^{\prime}(z)}$</code>, and the high-order iteration functions are merely a series in powers of <code>$u$</code>). 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).</p> http://mathoverflow.net/questions/37576/nth-order-generalizations-of-the-arithmetic-geometric-mean nth-order generalizations of the arithmetic-geometric mean J. M. 2010-09-03T04:01:48Z 2011-11-27T17:26:56Z <p>The arithmetic-geometric mean,</p> <p><code>$a_{k+1}=\frac{a_k+b_k}{2} \quad b_{k+1}=\sqrt{a_k b_k}$</code></p> <p>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$).</p> <p>I have been wondering if nth-order generalizations of the iteration,</p> <p><code>$a_{k+1}=\frac{a_k+b_k}{n} \quad b_{k+1}=\sqrt[n]{a_k b_k}$</code></p> <p>have ever been systematically studied. I've seen <a href="http://www.cecm.sfu.ca/personal/pborwein/PAPERS/P50.pdf" rel="nofollow">this paper</a> 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).</p> <p>Another possible generalization I thought involves $n$ variables and makes use of the <a href="http://en.wikipedia.org/wiki/Elementary_symmetric_polynomial" rel="nofollow">elementary symmetric polynomials</a>. To use $n=4$ as an example:</p> <p>$a_{k+1}=\frac{a_k+b_k+c_k+d_k}{4}$</p> <p>$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}}$</p> <p>$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}}$</p> <p>$d_{k+1}=\sqrt[4]{a_k b_k c_k d_k}$</p> <p>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?</p> <hr> <p><strong>EDIT</strong></p> <p>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</p> <p>$a_{k+1}=\frac{a_k+b_k+c_k+d_k}{4}$</p> <p>$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}}$</p> <p>$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}}$</p> <p>$d_{k+1}=\sqrt[4]{a_k b_k c_k d_k}$</p> <p>So, still the same question: is there a common limit, and if so, is the limit expressible in terms of known functions?</p> http://mathoverflow.net/questions/81553/orthogonal-polynomials-functions-on-the-interval-0-1-but-with-same-weight-as-ge/81574#81574 Answer by J. M. for Orthogonal polynomials/functions on the interval [0,1] but with same weight as Gegenbauer polynomials J. M. 2011-11-22T04:32:23Z 2011-11-22T13:12:06Z <p>There is <a href="http://www.cs.kuleuven.be/publicaties/rapporten/tw/TW534.pdf" rel="nofollow">this paper</a> and <a href="http://dx.doi.org/10.1016/j.cam.2011.10.006" rel="nofollow">this paper</a> 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.</p> <p>In the meantime, one can use the Stieltjes procedure to build up the recursion relations for these half range Gegenbauers. Letting</p> <p>$$\langle f(x),g(x) \rangle^{(\alpha)}=\int_0^1 (1-t^2)^{\alpha-1/2} f(t)g(t)\mathrm dt$$</p> <p>be the associated inner product, the Stieltjes procedure for generating monic orthogonal polynomials $\phi_k(x)$ uses the formulae</p> <p><code>\begin{align*}b_k&amp;=\frac{\langle x\phi_k(x),\phi_k(x)\rangle^{(\alpha)}}{\langle\phi_k(x),\phi_k(x)\rangle^{(\alpha)}}\\ c_k&amp;=\frac{\langle\phi_k(x),\phi_k(x)\rangle^{(\alpha)}}{\langle\phi_{k-1}(x),\phi_{k-1}(x)\rangle^{(\alpha)}}\end{align*}</code></p> <p>to give the coefficients $b_k,c_k$ for the recursion relation</p> <p>$$\phi_{k+1}(x)=(x-b_k)\phi_k(x)-c_k\phi_{k-1}(x)$$</p> <p>Here, the result</p> <p>$$\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)}$$</p> <p>is useful.</p> <hr> <p>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 <a href="http://math.stackexchange.com/questions/13174/34047#34047" rel="nofollow">here</a>, so I shall not repeat myself. Instead, I'll reproduce the <em>Mathematica</em> routine I gave in that answer:</p> <pre><code>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] ] </code></pre> <p>Here for instance is how to use <code>chebAlgo[]</code> to generate recursion coefficients for the monic half-range Chebyshev polynomials of the first kind:</p> <pre><code>With[{a = 0}, chebAlgo[Table[Gamma[(k + 1)/2] Gamma[a + 1/2]/Gamma[a + k/2 + 1], {k, 0, 10}]/2, Infinity]] // FullSimplify </code></pre> http://mathoverflow.net/questions/81576/sum-of-products-of-exponentials-and-polynomials/81583#81583 Answer by J. M. for Sum of products of exponentials and polynomials J. M. 2011-11-22T06:24:19Z 2011-11-22T06:24:19Z <p>You really can't do better than the polylogarithm and Lerch's transcendent as closed forms.</p> <p>For what it's worth, there are the <a href="http://functions.wolfram.com/ZetaFunctionsandPolylogarithms/PolyLog/03/01/01/02/0005/" rel="nofollow">explicit</a> <a href="http://functions.wolfram.com/ZetaFunctionsandPolylogarithms/LerchPhi/03/01/02/02/0001/" rel="nofollow">representations</a></p> <p><code>\begin{align*}\mathrm{Li}_{-n}(z)&amp;=\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)&amp;=a^n+\sum_{j=0}^n \binom{n}{j} \mathrm{Li}_{-j}(z) a^{n-j}\end{align*}</code></p> <p>all valid for $n$ a positive integer. Daunting, no?</p> <p>If you can't accept these closed forms, then you're better off staying with the sum representation that you have...</p> http://mathoverflow.net/questions/74478/on-a-polynomial-related-to-the-legendre-function-of-the-second-kind On a polynomial related to the Legendre function of the second kind J. M. 2011-09-04T02:27:04Z 2011-09-05T10:03:37Z <p>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.</p> <p>$Q_n(z)$ can be expressed in <a href="http://dlmf.nist.gov/14.7.E2" rel="nofollow">the following form</a>:</p> <p>$$Q_n(z)=P_n(z)\mathrm{artanh}\,z-W_{n-1}(z)$$</p> <p>where $W_{n-1}(z)$ can be expressed either <a href="http://dlmf.nist.gov/14.7.E4" rel="nofollow">as</a></p> <p>$$W_{n-1}(z)=\sum_{k=1}^n \frac{P_{k-1}(z) P_{n-k}(z)}{k}$$</p> <p>or <a href="http://dlmf.nist.gov/14.7.E3" rel="nofollow">as</a></p> <p>$$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$$</p> <p>where $H_k$ is the $k$-th harmonic number, $H_k=\sum\limits_{j=1}^k \frac1{j}$.</p> <p>My questions:</p> <ol> <li><p><em>Mathematica</em> returns a rather complicated expression for $W_{n-1}(z)$ involving the unknown solution of a certain recurrence (i.e. <code>DifferenceRoot[]</code>). Is there possibly a simpler form for this polynomial?</p></li> <li><p>Might there be a (hopefully simple) $n$-term recurrence that generates these polynomials?</p></li> </ol> <hr> <p><strong>Addendum</strong>:</p> <p>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.</p> <p>I now would like to expand my first question a bit: is it possible to express $W_n(z)$ as a <em>single</em> 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$?</p> <p>I would also like to consider an additional question: are the $W_n$ orthogonal polynomials with respect to some <em>weight function</em> $\omega(x)$ and an associated support interval $(a,b)$? That is, if</p> <p>$$\int_a^b\omega(t)W_j(t)W_k(t)\mathrm dt=0,\qquad j\neq k$$</p> <p>for some $\omega(x)$ and some interval $(a,b)$, what is this weight function and its support interval?</p> http://mathoverflow.net/questions/73234/roots-of-polynomial-with-matrix-coefficients/73327#73327 Answer by J. M. for roots of polynomial with matrix coefficients J. M. 2011-08-21T13:08:58Z 2011-08-21T13:08:58Z <p>Yes, there a number of numerical methods for finding the <em>solvents</em> of a polynomial with matrix coefficients. Dennis, Traub, and Weber in <a href="http://dx.doi.org/10.1137/0715034" rel="nofollow">this article</a> give some of the relevant theory, as well as some algorithms for finding so-called "dominant solvents" for your matrix polynomial (see also their <a href="http://dx.doi.org/10.1137/0713065" rel="nofollow">earlier article</a>). This <a href="http://dx.doi.org/10.1093/imanum/7.3.355" rel="nofollow">article</a> presents the use of Newton's method for solving matrix polynomials.</p> <p>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.</p> http://mathoverflow.net/questions/73027/how-many-curves-can-fit-on-a-sphere-without-intersecting How many curves can fit on a sphere without intersecting? J. M. 2011-08-17T03:07:14Z 2011-08-17T03:07:14Z <p>(This is more or less a generalized and more abstract version of <a href="http://math.stackexchange.com/questions/15801" rel="nofollow">this m.SE question</a>.)</p> <p>Let's recall two facts first:</p> <blockquote> <ol> <li><p>Any (arc-length parametrized) space curve is uniquely determined (up to rigid motions) by its curvature $\kappa(s)$ and its torsion $\tau(s)$.</p></li> <li><p>The curvature and torsion of a <em>unit</em>-spherical curve satisfies the equation $$\left(\tau\frac{\mathrm d}{\mathrm ds}\frac1{\kappa}\right)^2+\frac1{\kappa^2}=1$$</p></li> </ol> </blockquote> <p>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:</p> <blockquote> <p>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?</p> </blockquote> <p>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):</p> <p><code>\begin{align*}x&amp;=\cos\,\theta(t)\sin\,\varphi(t)\\y&amp;=\sin\,\theta(t)\sin\,\varphi(t)\\z&amp;=\cos\,\varphi(t)\end{align*}</code></p> <p>so we can also ask</p> <blockquote> <p>What conditions should $\theta(t)$ and $\varphi(t)$ satisfy so that $n$ copies of it can be placed on a sphere without intersecting?</p> </blockquote> <p>If there are better ways to characterize the conditions than these, I'd be interested in hearing about them.</p> <p>For convenience, let's consider only the unit sphere.</p> <p>One theorem(?) I have along these lines is</p> <blockquote> <p>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.</p> </blockquote> <p>but I haven't figured out how to generalize this.</p> http://mathoverflow.net/questions/72119/do-numerical-integration-with-fixed-abscissas/72148#72148 Answer by J. M. for do numerical integration with fixed abscissas J. M. 2011-08-05T05:34:46Z 2011-08-08T07:18:13Z <p>I have settled this question <a href="http://math.stackexchange.com/questions/55687/55742#55742" rel="nofollow">here</a>; in brief, a <a href="http://www.netlib.org/toms/655" rel="nofollow">FORTRAN implementation</a> of algorithms for solving this problem have been published in the Collected Algorithms of the ACM.</p> <hr> <p>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.</p> <p>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 <a href="http://www2.maths.ox.ac.uk/chebfun/" rel="nofollow"><code>chebfun</code> project</a> hinges on functions being nicely approximated by interpolating polynomials with abscissas that are transformed Chebyshev polynomial roots/extrema.)</p> <p>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...</p> http://mathoverflow.net/questions/72085/valid-to-use-all-wynn-extrapolated-values/72150#72150 Answer by J. M. for Valid to use all Wynn-extrapolated values? J. M. 2011-08-05T05:54:59Z 2011-08-05T05:54:59Z <p>There really isn't a theorem you can use, since most of the applicable theory is for sequences with <strong>known</strong> 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</p> <p>$$\varepsilon_{k+1}^{(n)}=\varepsilon_{k-1}^{(n+1)}+\frac1{\varepsilon_{k}^{(n+1)}-\varepsilon_{k}^{(n)}}$$</p> <p>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.</p> <p>For more practical advice, have a look at <a href="http://arxiv.org/abs/math/0306302" rel="nofollow">E.J. Weniger's survey paper</a>.</p> http://mathoverflow.net/questions/71905/numerical-beta-function/72149#72149 Answer by J. M. for Numerical Beta Function J. M. 2011-08-05T05:38:33Z 2011-08-05T05:38:33Z <p>There's nothing more straightforward than using the gamma function relationship, I believe (perhaps using <a href="http://my.fit.edu/~gabdo/gamma.txt" rel="nofollow">Lanczos's approximation</a> 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}$).</p> http://mathoverflow.net/questions/36207/finding-a-recursion-for-a-sum-of-legendre-polynomials Finding a recursion for a sum of Legendre polynomials J. M. 2010-08-20T17:45:42Z 2011-07-17T17:36:47Z <p>The polynomial</p> <p>$a_n(x):=P_n(x)-\frac{n-1}{n}P_{n-2}(x)$</p> <p>where $P_n(x)$ is a Legendre polynomial came up while I was investigating methods for estimating the error in Gaussian quadrature.</p> <p>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$...</p> <p>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?</p> http://mathoverflow.net/questions/35445/estimating-the-spectral-radius-of-a-matrix-noniteratively Estimating the spectral radius of a matrix, noniteratively J. M. 2010-08-13T06:31:10Z 2011-07-02T05:22:12Z <p>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:</p> <ol> <li>a similarity transformation to upper Hessenberg form (<code>$A=QHQ^T$</code>, <code>$Q$</code> orthogonal and <code>$H$</code> upper Hessenberg); and</li> <li>subtracting the identity multiplied by the mean of the eigenvalues from <code>$H$</code> (<code>$H'=H-\frac{trace(H)}{n}I$</code>, this corresponds geometrically to centering the eigenvalues around the origin of the complex plane).</li> </ol> <p>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 <code>$H'$</code> to its characteristic polynomial (equivalently, a similarity transformation of <code>$H'$</code> 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.</p> <p>My other naïve attempts at estimating the spectral radius include using <code>$||H'||_\infty$</code> 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.</p> <p>Is there a way to estimate the spectral radius more cheaply and noniteratively than actually computing eigenvalues?</p> http://mathoverflow.net/questions/39148/efficiently-computing-a-matrixs-induced-p-norm/39180#39180 Answer by J. M. for Efficiently computing a matrix's induced p-norm J. M. 2010-09-18T01:06:18Z 2011-05-13T04:10:04Z <p>Nicholas Higham gives an <a href="http://dx.doi.org/10.1007/BF01396242" rel="nofollow">algorithm</a> for <em>estimating</em> 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?</p> <hr> <p>(added 5/13/2011)</p> <p>I posted a <em>Mathematica</em> translation of Higham's <a href="http://people.sc.fsu.edu/~jburkardt/m_src/test_matrix/pnorm.m" rel="nofollow">original MATLAB code</a> <a href="http://math.stackexchange.com/questions/37913/37939#37939" rel="nofollow">here</a>.</p> http://mathoverflow.net/questions/64261/whats-the-difference-between-a-riemann-theta-and-a-siegel-theta-function What's the difference between a Riemann theta and a Siegel theta function? J. M. 2011-05-08T06:55:13Z 2011-05-08T20:11:41Z <p>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.</p> <p>A look at the <a href="http://dlmf.nist.gov/" rel="nofollow">DLMF</a> says that "the" multidimensional theta function is the <a href="http://dlmf.nist.gov/21.2" rel="nofollow">Riemann theta function</a>,</p> <p>$$\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)}}$$</p> <p>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 <a href="http://www.ams.org/journals/mcom/2004-73-247/S0025-5718-03-01609-0/S0025-5718-03-01609-0.pdf" rel="nofollow">nice implementation</a> for numerically computing this function, but it is only available in Maple and Java. Since I work in <em>Mathematica</em>, I decided to see if there was already something in <em>Mathematica</em> before going through the trouble of translating code.</p> <p>Peering at the docs for <em>Mathematica</em> netted the <a href="http://reference.wolfram.com/mathematica/ref/SiegelTheta.html" rel="nofollow">Siegel theta function</a>. However, looking at how the Siegel theta function was defined, this and the Riemann theta function <em>seem</em> to be the same thing!</p> <p>Not having Maple to check if the results from their implementation of Riemann theta and the results from <em>Mathematica</em>'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?</p> http://mathoverflow.net/questions/64302/numerical-differentiation-what-is-the-best-method/64306#64306 Answer by J. M. for Numerical Differentiation. What is the best method? J. M. 2011-05-08T15:33:31Z 2011-05-08T15:33:31Z <p>If your function is badly behaved (e.g. noisy, very oscillatory), no method will perform properly (differentiation is numerically very <em>unstable</em>). That being said, for "nice functions", I have good experience with polynomial (Richardson) extrapolation methods. <a href="http://dx.doi.org/10.1007/BF02166671" rel="nofollow">This paper</a> and <a href="http://dx.doi.org/10.1007/BF01933664" rel="nofollow">this paper</a> give hints on how you might write your own implementation. I will note that this is the method implemented in the <a href="http://nag.co.uk/numeric/numerical_libraries.asp" rel="nofollow">NAG numerical libraries</a> (with of course a few wrinkles of their own).</p> <p>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:</p> <p>$$f^\prime(x)=\frac1{2\pi i}\oint_\gamma \frac{f(t)}{(t-x)^2}\mathrm dt$$</p> <p>where it is up to you to choose a suitable counterclockwise contour $\gamma$ (a circle is customary); the other is to use the <a href="http://www.jstor.org/pss/2589707" rel="nofollow">"Lanczos derivative"</a>:</p> <p>$$f^\prime(x)=\lim_{h\to 0}\frac{3}{2h^3}\int_{-h}^h t\;f(x+t)\mathrm dt$$</p> <p>where you either will have to experiment with an appropriate step size $h$, or use some extrapolative procedure.</p> <p>You will have to experiment with your computing environment to choose.</p> http://mathoverflow.net/questions/50544/some-dirichlet-series-questions Some Dirichlet series questions. J. M. 2010-12-28T06:31:24Z 2011-01-13T00:22:13Z <p>I asked <a href="http://math.stackexchange.com/questions/15350" rel="nofollow">this question</a> on m.SE in an attempt to find out the right words to say for these questions I am about to ask.</p> <p>In his great answer, Matthew Emerton explained that (cuspidal) automorphic <em>L</em>-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:</p> <p>1) Are there Dirichlet series that cannot be classified as (cuspidal) automorphic <em>L</em>-functions, yet still possess a critical line of nontrivial zeroes?</p> <p>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</p> <p>$$g(\sigma+it)=z(t)\exp(-i\vartheta(t))$$</p> <p>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</p> <p>2) Do all (cuspidal) automorphic <em>L</em>-functions have a "Riemann-Siegel" decomposition? If not, what restrictions are there for them to possess such a decomposition?</p> <p>My motivation is more of curiosity than anything else. Hopefully this is not too elementary a question!</p> http://mathoverflow.net/questions/50343/what-would-you-want-to-see-at-the-museum-of-mathematics/50346#50346 Answer by J. M. for What would you want to see at the Museum of Mathematics? J. M. 2010-12-25T16:28:51Z 2010-12-25T16:28:51Z <p>Sculptures of surfaces would be lovely.</p> http://mathoverflow.net/questions/40876/second-order-taylor-expansion-to-solve-system-of-equations/50168#50168 Answer by J. M. for Second order Taylor expansion to solve system of equations J. M. 2010-12-22T15:56:43Z 2010-12-22T15:56:43Z <p>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 <a href="http://dx.doi.org/10.1145/3147.3162" rel="nofollow">multivariate version of the Halley iteration</a> has already been developed decades ago. Maybe this might be of use.</p> http://mathoverflow.net/questions/47658/why-not-evaluate-integrals-using-ode-solvers/47776#47776 Answer by J. M. for Why not evaluate integrals using ODE-solvers? J. M. 2010-11-30T09:42:14Z 2010-11-30T10:05:53Z <p>Here's my take on the matter: the difference of philosophy between quadrature routines and ODE solving routines, I believe, is this:</p> <blockquote> <p>Extrapolation is riskier than interpolation.</p> </blockquote> <p>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.</p> <p>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 <em>extrapolation is inherently unstable</em>. 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.</p> <p>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.</p> <p>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 <em>in a relatively tiny interval within the interval of integration</em>, 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.</p> http://mathoverflow.net/questions/46553/fast-trace-of-inverse-of-a-square-matrix/46620#46620 Answer by J. M. for Fast trace of inverse of a square matrix J. M. 2010-11-19T12:35:07Z 2010-11-19T12:35:07Z <p>Given that the poster has specified that his matrix is symmetric, I offer a general solution and a special case:</p> <ol> <li><p>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.</p></li> <li><p>If the matrix is positive definite as well, first perform a Cholesky decomposition. Then there are <a href="http://books.google.com/books?id=myzIPBwyBbcC&amp;pg=PA119" rel="nofollow">methods</a> for generating the diagonal elements of the inverse.</p></li> </ol> http://mathoverflow.net/questions/44326/most-memorable-titles/44330#44330 Answer by J. M. for Most memorable titles J. M. 2010-10-31T14:46:55Z 2010-10-31T14:46:55Z <p>The AKS paper <a href="http://www.cse.iitk.ac.in/users/manindra/algebra/primality_v6.pdf" rel="nofollow">Primes is in P</a> is a pretty memorable title for me.</p> http://mathoverflow.net/questions/43943/square-root-of-non-positive-definite-matrix/43955#43955 Answer by J. M. for Square root of non-positive definite matrix J. M. 2010-10-28T09:08:21Z 2010-10-28T09:08:21Z <p>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$)</p> <p>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!</p> <p>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).</p> <p>I have only given a brief sketch, since what I really should be doing is to point out this <a href="http://eprints.ma.man.ac.uk/1193/01/covered/MIMS_ep2008_56.pdf" rel="nofollow">paper</a> by Nick Higham, which discusses the nuances of Cholesky decomposition with symmetric pivoting for symmetric positive semidefinite matrices.</p> http://mathoverflow.net/questions/43514/how-do-eigenvectors-and-eigenvalues-change-when-we-remove-a-row-column-pair-of-a/43516#43516 Answer by J. M. for How do eigenvectors and eigenvalues change when we remove a row/column pair of a matrix? J. M. 2010-10-25T13:36:05Z 2010-10-25T13:36:05Z <p>The word you're looking for is <em>downdating</em>, and I cannot do better than to point out <a href="http://onlinelibrary.wiley.com/doi/10.1111/1467-9574.00174/abstract" rel="nofollow">these</a> <a href="http://www.computer.org/portal/web/csdl/doi/10.1109/MCSE.2006.36" rel="nofollow">two</a> survey papers, and this <a href="http://dx.doi.org/10.1016/S0262-8856%2802%2900114-2" rel="nofollow">article</a>. 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}$.</p> http://mathoverflow.net/questions/132186/regarding-a-paper-by-paul-a-clement-on-tridiagonal-matrices Comment by J. M. J. M. 2013-05-30T19:04:25Z 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 &quot;better&quot;...) http://mathoverflow.net/questions/132186/regarding-a-paper-by-paul-a-clement-on-tridiagonal-matrices/132261#132261 Comment by J. M. J. M. 2013-05-30T19:03:06Z 2013-05-30T19:03:06Z As you have linked to the paper by Taussky and Todd, you might want to mention <a href="http://www-math.mit.edu/~edelman/homepage/papers/kac.pdf" rel="nofollow">www-math.mit.edu/~edelman/homepage/papers/kac.pdf</a> as well... http://mathoverflow.net/questions/32099/what-is-lagrange-inversion-good-for/84131#84131 Comment by J. M. J. M. 2011-12-23T01:50:10Z 2011-12-23T01:50:10Z That solution gives a complex root, not a real one. http://mathoverflow.net/questions/84079/how-to-find-the-region-of-convergence-of-this-series-using-the-theory-of-continue/84082#84082 Comment by J. M. J. M. 2011-12-23T01:36:56Z 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. http://mathoverflow.net/questions/84108/the-riemanns-zeta-function-represented-as-a-continued-fraction-and-a-question-of Comment by J. M. J. M. 2011-12-23T01:17:54Z 2011-12-23T01:17:54Z @Noam: it's somewhat standard in CF literature. I'm told it's originally Gauss's. http://mathoverflow.net/questions/36272/when-is-arctan-a-rational-multiple-of-pi Comment by J. M. J. M. 2011-12-18T10:47:19Z 2011-12-18T10:47:19Z The article @Doug linked to is also at <a href="http://www.oberlin.edu/faculty/jcalcut/gausspi.pdf" rel="nofollow">oberlin.edu/faculty/jcalcut/gausspi.pdf</a> ; see also <a href="http://www.oberlin.edu/faculty/jcalcut/arctan.pdf" rel="nofollow">oberlin.edu/faculty/jcalcut/arctan.pdf</a> . http://mathoverflow.net/questions/3440/spiral-of-theodorus/3444#3444 Comment by J. M. J. M. 2011-12-16T14:32:34Z 2011-12-16T14:32:34Z <a href="http://www.uni-graz.at/~gronau/monthly230-237.pdf" rel="nofollow">uni-graz.at/~gronau/monthly230-237.pdf</a> is a newer paper discussing the analytic continuation of the Theodorus spiral. http://mathoverflow.net/questions/83608/symbolic-diagonalization-of-a-matrix Comment by J. M. J. M. 2011-12-16T12:31:42Z 2011-12-16T12:31:42Z <i>Mathematica</i> and Maple sidestep the problem of &quot;symbolic expression&quot; 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 <code>$\geq 5$</code>, you need theta functions. Which is a pretty deep rabbit hole... http://mathoverflow.net/questions/40802/computation-of-a-drazin-inverse Comment by J. M. J. M. 2011-12-15T01:51:46Z 2011-12-15T01:51:46Z Eep, sorry for not checking that first DOI! Here's the first one: <a href="http://www.springerlink.com/content/u528652t80v96515/" rel="nofollow">springerlink.com/content/u528652t80v96515</a> http://mathoverflow.net/questions/40802/computation-of-a-drazin-inverse Comment by J. M. J. M. 2011-12-14T13:28:33Z 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 <a href="http://dx.doi.org/10.1007/BFb0120751" rel="nofollow">dx.doi.org/10.1007/BFb0120751</a> and <a href="http://dx.doi.org/10.1137/0131057" rel="nofollow">dx.doi.org/10.1137/0131057</a> might be applicable to your problem. http://mathoverflow.net/questions/22811/upper-bound-of-period-length-of-continued-fraction-representation-of-very-composi/23014#23014 Comment by J. M. J. M. 2011-12-10T11:33:17Z 2011-12-10T11:33:17Z For cycle detection, either of Floyd's or Brent's algorithms should be suitable. http://mathoverflow.net/questions/83023/how-to-do-integrals-involving-two-bessel-functions-and-another-function Comment by J. M. J. M. 2011-12-10T04:26:10Z 2011-12-10T04:26:10Z ren, you might want to also try asking this question on <a href="http://math.stackexchange.com" rel="nofollow">math.stackexchange.com</a> ; 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. http://mathoverflow.net/questions/72616/solving-for-moore-penrose-pseudo-inverse/72627#72627 Comment by J. M. J. M. 2011-12-09T04:22:36Z 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&amp;1\\1&amp;0\end{pmatrix}$. http://mathoverflow.net/questions/47658/why-not-evaluate-integrals-using-ode-solvers/47662#47662 Comment by J. M. J. M. 2011-12-09T03:57:14Z 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 <a href="http://dx.doi.org/10.1007/BF02165405" rel="nofollow">dx.doi.org/10.1007/BF02165405</a> and <a href="http://dx.doi.org/10.1007/BF02167558" rel="nofollow">dx.doi.org/10.1007/BF02167558</a> for instance. There is now the algorithm due to Carlson: <a href="http://dx.doi.org/10.1007/BF02198293" rel="nofollow">dx.doi.org/10.1007/BF02198293</a> , which is slightly less complicated, but I haven't done a comparison of this versus plain quadrature with, say, Gauss-Kronrod. http://mathoverflow.net/questions/82739/approximation-algorithmic-considerations/82767#82767 Comment by J. M. J. M. 2011-12-09T03:43:40Z 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&#233; 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.