15

In a single application, it is rather likely that you only will need a small subset of all the possible extremes of the generalized hypergeometric function. It is a very general function, after all. Having an idea about the range of $z$ and of the parameters $a_i, b_i$ would allow giving more specific advice. In general, the standard method, assuming $p \le ...


12

I've written my own integrator, quadcc, which copes substantially better than the Matlab integrators with singularities, and provides a more reliable error estimate. To use it for your problem, I did the following: >> lambda = 0.00313; kappa = 0.00825; nu = 0.33; >> x = 10; >> E = @(r) r.^4.*(lambda*sqrt(kappa^2 + r.^2)).^(-nu-5/2) .* ...


11

As Pedro points out, Levin-type methods are the best established methods for these kinds of problems. Do you have access to Mathematica? For this problem, Mathematica will detect and use them by default: In[1]:= e[r_] := r^4 (l Sqrt[k^2 + r^2])^(-v - 5/2) BesselK[-v - 5/2, l Sqrt[k^2 + r^2]] In[2]:= {l, k, v} = {0.00313, 0.00825, 0.33}; In[3]:= Block[{...


11

This is a root-finding problem for an analytic function over a single dimension. One standard technique is to approximate $F(k)$ as a polynomial $F(k) \approx c_0 + c_1 k + c_2 k^2 + \cdots$ using a Chebyshev approximation, and to compute the roots of the polynomial semi-analytically, e.g. by setting up the companion matrix and computing its eigenvalues. (...


10

Doing one-off best rational approximations can often be accomplished by "manual" iterations of the Remez algorithm: interpolate a rational approximation with (relative or absolute) alternating sign errors at an initial guess for interpolation points, locate one (or more) points where the actual error exceeds that of the interpolation points and pivot (...


9

Do a polynomial expansion (a la l'Hopital's rule) for both enumerator and denominator and you get a rational function that, for small $x$, will approximate the function well. As an example: $$ \frac{\sin ax}{\sin x} \approx \frac{ax-\tfrac 1{3!}(ax)^3 + \ldots}{x-\tfrac 1{3!}x^3 + \ldots} = \frac{a-\tfrac 1{3!}a^3x^2 + \ldots}{1-\tfrac 1{3!}x^2 + \...


9

The integral in question is also known as the Boys function, after the British chemist Samuel Francis Boys who introduced its use in the early 1950s. A few years ago, I needed to compute this function in double precision, as fast as possible but accurately. I managed to achieve a relative error on the order of $10^{-15}$ across the entire input domain. It ...


9

There is a prime missing both in your code and in the expression in your question. In the original paper, the expression is: $$ \log(x/x_0) + (x-x_0) \frac{\sum^\prime_{0\leq j\leq n}p_j T_j^*(x/6)}{\sum^\prime_{0\leq j\leq n}q_j T_j^*(x/6)}. $$ This primed sum is very common and standard when dealing with Chebyshev series: it means that the first term of ...


8

Cube roots are not nearly as important as square roots (e.g., for normalizing vectors), so that might be why they are discussed much less. In general, if you apply Newton's method to $x^\alpha-\beta$, you get the iteration $$ (1-\alpha^{-1})x+\frac\beta\alpha x^{1-\alpha}, $$ so you just need to pick the equation so that $\alpha$ is negative integer, it's ...


8

I will complement @Richard Zhang 's answer (+1) with a python implementation of his suggested approach. The MATLAB package Chebfun has been partially ported in python. Actually there are two versions available: chebpy and pychebfun. Here's an implementation of the root finding procedure with pychebfun (the approach is similar with chebpy) from scipy....


7

You can use the following recurrence relation $\mathop{S}\nolimits\left(n,k\right)=k\mathop{S}\nolimits\left(n-1,k\right)+\mathop{S}\nolimits\!\left(n-1,k-1\right)$ see DLMF equation 26.8.E22 and build the (triangular) table.


7

I've encountered a similar problem before, but I was more worried about speed than accuracy (see paper here). If your angle $\vartheta$ is the result of an $\arccos(\cdot)$, which is often the case in geometric computations, you can use Chebyshev polynomials which are defined as $$T_k(x) = \cos(k\arccos(x)), \quad \mbox{or} \quad T_k(\cos(\vartheta)) = \cos(...


7

The problem is that np.cos(t) and np.sqrt(t) generate arrays with the length of t, whereas the second row ([0,1]) maintains the same size. To use np.vectorize with your function, you have to define the output type, and np.vectorize isn't really meant as a decorator except for the simplest cases. In this way however you can generate the function with the ...


6

What have you tried so far? This completely naive implementation manages to compute 7 (maybe 7.5) digits in 2.5 seconds on my laptop: #include <iostream> #include <complex> #include <cmath> #include <iomanip> int main () { const double alpha = 1; std::cout.precision(16); std::complex<double> sum = 0; for (unsigned int ...


6

There's a GPL'd C library, ANANT - Algorithms in Analytic Number Theory by Linas Vepstas, which includes multiprecision implementation of the polylogarithm, building on GMP. From its README file: This project contains ad-hoc implementations of assorted analytic functions of interest in number theory, including the gamma function, the Riemann zeta ...


6

For the Hankel transform, one can classify the methods into four major groups: Numerical quadrature-based. Fourier-based ones. Asymptotic expansion of Bessel into sines and cosines. Projection-slice methods. The following paper gives a nice overview of these methods (types of methods): M. J. Cree and P. J. Bones, "Algorithms to numerically evaluate the ...


6

You're probably better served writing a C wrapper for the Fortran implementation you linked to: Colavecchia, F. D., Gasaneo, G., "f1: a code to compute Appell's F1 hypergeometric function", Computer Physics Communications, Volume 157, Issue 1, p. 32-38 (2004), found at http://cpc.cs.qub.ac.uk/summaries/ADSJ. The R package appell wraps that implementation of ...


6

The appropriate and fastest library depends on several things. Which Bessel functions (only J, Y & Hankel or modified Bessel functions I & K too), for which types of arguments (real or complex, integer, fractional or general order)? Amos's libraries are written in Fortran-77 (there are Fortran-90 coverted versions of TOMS 644 on a mirror of Alan ...


6

You can avoid overflow by rewriting the sum appropriately. Say you have a sum $$ X = \sum_{i=1}^{n} x_i, $$ where some $x_i$'s are positive and large enough to cause overflow, but $\log x_i$ are of reasonable magnitude. If you find the greatest $\log x_M \geq \log x_i$, then you can rewrite the sum as $$ \log X = \log \sum_i e^{\log x_i} = \log x_M + \log \...


6

After some searching I found two formulas to solve the accuracy problem. The two modifications are In range $|k| \ge 40$ the (asymptotic) series from http://functions.wolfram.com/08.01.06.0010.02 is used (it is valid for $|k|>1,\,$ but convergence is suboptimal for small $k$): \begin{multline*} E(z)=\sqrt{-z} +\frac{\ln(-z)}{\sqrt{-z}}\sum_{n=0}^\infty \...


6

First of all, from the first paragraph of your attempts at a solution, I assume that the $z_j$ are non-negative? In that case, the integrand has no real problematic points (it's monotonous, decreasing). The only difficulty is the infinite integration range. And then it all depends on the accuracy you want... Do you want machine epsilon precision (or there ...


5

If you don't have access to Mathematica, you could write a Levin-type (or other specialized oscillatory) method in Matlab as Pedro suggests. Do you use the chebfun library for Matlab? I just learned it contains an implementation of a basic Levin-type method here. The implementation is written by Olver (one of the experts in the oscillatory quadrature field)....


5

Using the half angle formulas you can convert the exponent into the form $$a + b \cos 2\theta + c \sin 2\theta$$ which then integrates nicely into Bessel functions. Mathematica gives $$\int_0^{2\pi}e^{a + b \cos 2\theta + c \sin 2\theta} d\theta = e^a \pi I_0\left(\sqrt{b^2 + c^2}\right)$$


5

You might take a look at Numerical Methods for Special Functions by Amparo Gil, Javier Segura, and Nico M. Temme.


5

As you can see in this plot of $\log|f(\lambda)|$, $$ f(\lambda) = J_{\lambda-1}(1) - 2J_\lambda(1) -J_{\lambda+1}(1), $$ the roots $\lambda_k$ are really regular, and are approximately equal to $-k$ (starting from $k\geq0$, $\lambda_0=1.23219$ is an exception). So the way to get the $k$-th eigenvalue ($k\geq1$), is to bracket the root in $[-k-\frac12,-k+\...


5

Looking at what the python scipy library does for its special functions, the polygamma is found by returning the digamma if the zeroth derivative is requested, otherwise return $(-1)^{n+1}\Gamma(n+1)\zeta(n+1,z)$ where $\Gamma$ is the gamma function and $\zeta$ the two argument Riemann zeta function. Assuming that this identity holds for the complex numbers, ...


5

I hope I understood the question correctly. They try to compute exactly the same thing, so they really are equivalent. I'll use Chebyshev polynomials because they are easy to analyze. Given a function $f(x)$ on $[-1,1]$, the spectral interpolant is the truncation of $$ \begin{aligned} f(x) &= \sum_{n\geq0} \bar a_n T_n(x), \\ \bar a_n &= \frac{1+[n&...


5

Step 1: analyze the recurrence I implemented the recurrence in Python, using basic numpy and scipy.special for the erf function. Why? Because it is simple and cheap (for reasonable $p$). import numpy as np import scipy.special as sp def F1(z): return np.sqrt(np.pi/2)*(1.0+sp.erf(z/np.sqrt(2))) def F2(z): return F1(z)*z + np.exp(-z**2/2.0) def F(...


5

Complementing the answers, I would like to say that directly discretizing the differential equation might work really well. I used that approach in a previous question. I used Finite Differences and used an eigensolver to find the roots. That, besides a perturbation analysis followed by Newton's method. The suggestion there was to use Chebyshev ...


Only top voted, non community-wiki answers of a minimum length are eligible