Logo

Numerical examples




A dynamical computation of a multiple integral involved in the neutron star theory

This study consists in computing dynamically using the CADNA library an integral involved in the neutron star theory. Massive stars collapse with an explosion ejecting all the matter, except a central residue: a neutron star or a black hole, depending on the details of the scenario. Neutron stars are mainly done of neutrons with a typical proportion of protons smaller than 10%. With time, equilibrium between protons and neutrons takes place. Time to reach this equilibrium is obtained by calculating high order integrals involving several variables for any particle appearing in the reactions.


The integral, which has been provided by Loïc Villain (LUTH, Observatoire de Paris-Meudon, France), is the following:

\begin{displaymath}I \left( \varepsilon ,v \right) = \frac{1}{g \left( \varepsil... ...^{ \infty }_0 dp\ H \left( n, p, \theta, \varepsilon, v \right)\end{displaymath}

where $(\varepsilon,v) \in [10^{-4},10^{4}]\times[10^{-4},10^{3}]$ and $g$ is a normalization function.

$H \left( n, p, \theta, \varepsilon, v \right) = f(z)\Gamma(n-\varepsilon-z)+f(-z)\Gamma(n-\varepsilon+z) -f(z)\Gamma(n+\varepsilon-z)-f(z)\Gamma(n+\varepsilon+z)$,

with $z = \sqrt{p^2 + y\left( v, \theta \right) ^2},\ f(x)=\frac{1}{e^{x}+1}, \ \Gamma(x) = \frac{x}{e^{x}-1} \ $ and $\ y(v, \theta )=v \sin ( \theta )$.


To compute $I \left( \varepsilon ,v \right)$, several numerical problems have to be faced:

  • The integral has two infinite bounds. Each integral of the form $\int_0^\infty$ is replaced by $\sum_{j=0}^k \int_{jL}^{(j+1)L}$. One problem is the choice of an optimal value for $k$.
  • $\Gamma(x) = \frac{x}{e^{x}-1}$ generates cancellations if $x\approx 0$.
    Therefore a series expansion of $\Gamma(x)$ is used: $\Gamma(x) \approx \frac{1}{1+ \frac{x}{2} + ... +\frac{x^{n-1}}{n!}}$ One problem is the choice of an optimal value for $n$.
  • $I \left( \varepsilon ,v \right)$ is the composition of three one-dimensional integrals. A quadrature method, such as Gauss-Legendre method, is used to compute each one-dimensional integral. The truncation error generated by this quadrature method is another problem.


$I \left( \varepsilon ,v \right)$ is expressed as one-dimensional integrals:

\begin{displaymath}I\left( \varepsilon ,v \right) = \frac{1}{g \left( \varepsil... ...^{ \infty }_0 H \left( n, p, \theta, \varepsilon, v \right) dp.\end{displaymath}


\begin{displaymath}I\left(\varepsilon ,v \right) = \frac{1}{g \left( \varepsilon... ...nt^{ \frac{ \pi }{ 2 } }_0 f_1(\theta , \varepsilon, v) d\theta\end{displaymath}


\begin{displaymath}{\rm with} \ \ \displaystyle f_1(\theta , \varepsilon, v) = sin(\theta) \int^{ \infty }_0 f_2(n,\theta , \varepsilon, v) dn\end{displaymath}


\begin{displaymath}{\rm and} \ \ \displaystyle f_2(n,\theta, \varepsilon, v) =n... ...nt^{ \infty }_0 H \left(n, p, \theta, \varepsilon, v \right) dp\end{displaymath}

Each one-dimensional integral is evaluated by two iterative processes:

  • A sequence $(B_n)$ of integrals on a finite interval with an increasing upper bound:

    \begin{eqnarray*} B_n & = & \int_0^{nL} f(t)dt\\ B_{n+1} & = & B_{n} + \underbrace{\int_{nL}^{(n+1)L}f(t)dt}_{a_n} \end{eqnarray*}

    The sequence $(B_n)$ is computed until the difference $B_n - B_{n+1}$ is not significant.
  • Each integral $a_n$ is computed using Gauss-Legendre method. $[nL,(n+1)L]$ is partitioned in $2^p$ sub-intervals of length $\frac{L}{2^p}$. The integral on each sub-interval is approximated using the classical Gauss-Legendre method. Let $A_n^p$ be the approximation of $a_n$ using Gauss-Legendre method on $2^p$ sub-intervals of $[nL,(n+1)L]$. The sequence $A_n^p$ $(p=1,2,...)$ is computed until the difference $A_n^p - A_n^{p+1}$ is not significant.

Using this strategy, the exact singificant digits of the approximation of $I \left( \varepsilon ,v \right)$ obtained are those of the exact value of the integral, up to 6.


To get the FORTRAN source code.

More information can be requested to the CADNA team
Thanks to Baptiste Mary for the CADNA logo