2
$\begingroup$

I am working on a math problem as part of my Physics project, and I need to numerically evaluate the integral

$\int_0^{\infty}\frac{f(x)}{(1+ax)^n}dx$

where $a$ is some positive constant and $1 < n < 2$, also $f(x)$ can be assumed to be slowly varying.

The problem is that I don't actually have an analytical expression for $f(x)$, I only have samples of $f(x)$ for a finite range of $x$. And assuming that $f(x)$ is slowly varying, the fractional error in stopping the integration at some limit $x'$ is $\sim (1+ax')^{1-n}$. This implies that for a given allowed tolerance, $x'$ becomes really large as $n$ approaches 1. This is a problem as sampling $f(x)$ for large $x$ costs a lot of computation time.

What are some ways I can reduce the computation time needed without losing accuracy?

$\endgroup$
14
  • $\begingroup$ Do you know the value of $f(x)$ for all $x\in(0,x')$? That is, are you given a continuous set of points to work with? $\endgroup$
    – Simply Beautiful Art
    Feb 9, 2017 at 21:50
  • $\begingroup$ @SimplyBeautifulArt I can interpolate $f(x)$ from a finite set of points for $x<x'$. $\endgroup$
    – David Young
    Feb 9, 2017 at 21:53
  • $\begingroup$ @Dr.MV: you are right, but I think there is a numerical work-around for that, please see my answer below. $\endgroup$
    – Jack D'Aurizio
    Feb 9, 2017 at 22:44
  • $\begingroup$ It would be useful to have some authority on numerical analysis here. Is there a way to summon sir Claude Leibovici, by your knowledge? $\endgroup$
    – Jack D'Aurizio
    Feb 9, 2017 at 22:45
  • 1
    $\begingroup$ You might want to look into the Generalized Gauss-Laguerre-quadrature. $\endgroup$
    – flawr
    Feb 9, 2017 at 22:46

1 Answer 1

Reset to default
0
$\begingroup$

The definition of the incomplete Gamma function is $$ \Gamma(s,z_0) = \int_{z_0}^{\infty} z^{s-1} e^{-z} \, dz $$ (a complex contour integral if $z_0$ is complex). Setting $s=1-n$ and $z_0=-ib/a$ and setting the integration path to $z(t) = -ib/a - ibt$ (for $t$ from 0 to $\infty$) results in $$ \Gamma(1-n, -ib/a) = -\frac{a^n}{b^{n-1}} i^{n+1} \exp(ib/a) \int_0^{\infty} \frac{\exp(ibx)}{(1 + ax)^n} \, dx $$ or $$ \int_0^{\infty} \frac{\exp(ibx)}{(1 + ax)^n} \, dx = -\frac{b^{n-1}}{a^n} (-i)^{n+1} \exp(-ib/a) \Gamma(1-n, -ib/a). $$ Use $(-i)^{n+1} = \exp((n+1) \log(-i))$ with branch cut of the negative real axis to evaluate $(-i)^{n+1}.$ With standard branch cut of the negative real axis, $\log (-i) = -i\pi/2.$

So if you can express your function $f$ as a linear combination of sinusoids, you can use the incomplete Gamma function to evaluate this integral. Evaluating the incomplete Gamma function is somewhat challenging, but most numerical analysis software will do it.

$\endgroup$
4
  • $\begingroup$ Thanks for the interesting suggestion and sorry for not having responded. I think this would be useful if I have $f(x)$ as a Fourier series then I can integrate each term, but I am not sure how this can be applied to my problem where I only have $f(x)$ for a finite domain. $\endgroup$
    – David Young
    Feb 16, 2017 at 17:20
  • $\begingroup$ You would have to do a non-linear least squares fit (e.g., Levenberg-Marquardt) to your data over the finite interval with a linear combination of sinusoids (meaning $f(x)=\sum_{i=1}^N a_i \sin(2\pi b_i x + c_i)$, where the $a_i,$, $b_i,$ and $c_i$ are the unknowns) to get the best fit to your data. And then hope that it extrapolates well to the data beyond the finite interval. The best way to ensure that it extrapolates well is to use the smallest $N$ possible (to avoid "overfitting") while fitting the data in the finite interval well. This might be a topic for another question. $\endgroup$
    – J. Heller
    Feb 16, 2017 at 17:53
  • $\begingroup$ I am not sure if it would extrapolate well beyond the interval. I have worked out that a specific case of my problem can be modelled roughly with the integral $\int_0^{x_{lim}}\frac{sin(x)}{(1+400x)^{1.8}}dx$, where $x_{lim} \sim 1$. $\endgroup$
    – David Young
    Feb 16, 2017 at 18:25
  • $\begingroup$ BTW, the Wikipedia article on Fourier series gives the Fourier series approximating a function on an arbitrary interval $[x_0,x_0+P]$. You could use the Fourier series to get a good initial guess for non-linear least squares fitting (using the Fourier series terms with largest magnitudes). $\endgroup$
    – J. Heller
    Feb 16, 2017 at 18:33

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .