Definition of $\Gamma (a,x)$

The big_gamma function, $\Gamma (a,x)$, is given in terms of a continued fraction expansion:

$$\Gamma(a, x) e^{-x}x^a \frac{1} {x + \frac{1 - a} {1 + \frac{1} {x + \frac{2-a... ... + \frac{2} {x + \frac{3-a} {1 + \frac{3} {x + \frac{4 - a} {1 + \ldots}}}}}}}}$$ =

$$e^{-x}x^a \left(\frac{1}{x+}\, \frac{1-a}{1+} \, \frac{1}{x+} \, ... ...-a}{1+} \frac{2}{x+}\, \cdots \, \frac{n-a}{1+} \, \frac{n}{x+}\, \cdots\right)$$ = (2.27)

Continued fraction expansions may look dreadful but they converge extremely rapidly compared to power series expansions. For this reason they've been of great interest to numerical analysts and, of course, to computational scientists too. In case of the incomplete gamma function the continued fraction expansion converges best where the series expansion does not, so there is also a very nice complementarity here.

Unlike a power series expansion the continued fraction expansion seemingly cannot be evaluated in steps from ``left to right'', i.e., from the lowest terms, gradually adding higher and higher terms to the sum. Here it looks like we have to start from the end, sic! That is, we have to guess the last term, and then gradually build up the whole fraction starting from the last term in the expansion and moving towards the beginning.

But it can be done the right way, i.e., from ``left to right''. Consider the following generic continued fraction expansions:

$$b_0 + \frac{a_1}{b_1 +}\,\frac{a_2}{b_2 +}\, \cdots \, \frac{a_{n-2}}{b_{n-2}+}\, \frac{a_{n-1}}{b_{n-1}}$$ f_n-1(x) =

$$\frac{A_{n-1}}{B_{n-1}}$$ =

$$b_0 + \frac{a_1}{b_1 +}\,\frac{a_2}{b_2 +}\, \cdots \, \frac{a_{n-1}}{b_{n-1}+}\, \frac{a_{n}}{b_{n}}$$ f_n(x) =

$$\frac{A_n}{B_n}$$ = (2.28)

The question is, how A_n and B_n relate to A_n-1, B_n-1, A_n-2, B_n-2, etc.

We can help ourselves with Maxima   in this task. First we define f₁ and extract its numerator and denominator:

Now we'll do likewise with f₂:

Observe that without asking Maxima has grouped biga[2] in a way that shows the presence of biga[1] in it! Let us now evaluate f₃:

Let us list our obtained expressions for A_i and B_i in order, and attempt some substitutions:

In order to carry out substitutions we use Maxima  function ratsubst, which takes three arguments: the name of a new variable, its value, and the expression into which the variable should be substituted:

What we seem to be getting is the following relation:

A₃ = b₃ A₂ + a₃ A₁ (2.29)

Now let us have a closer look at B_i:

Here we got that

B₃ = b₃ B₂ + a₃ B₁ (2.30)

Could it be that the general formula is:

  

A_n = b_n A_n-1 + a_n A_n-2 (2.31)

B_n = b_n B_n-1 + a_n B_n-2 (2.32)

We can check it by evaluating f₄:

and performing appropriate substitutions:

Well, formulas 2.31 and 2.32 clearly hold. But the recurrence relationships for A_n and B_n are identical, so clearly, if they are not to be the same their start-up values must be different. What are they? Let us have a look at A₁ first:

Matching this against equation 2.31 yields:

A₁ = b₁ A₀ + a₁ A_-1 = b₁ b₀ + a₁ 1

hence

$$A_0 = b_0\qquad\textrm{and}\qquad A_{-1} = 1$$ (2.33)

Similarly for B₁:

Matching this against equation 2.32 yields:

B₁ = b₁ B₀ + a₁ B_-1 = b₁ 1 + a₁ 0

hence

$$B_0 = 1\qquad\textrm{and}\qquad B_{-1} = 0$$ (2.34)

The evaluation of a continued fraction expansion now looks much easier. All that we have to do is to

1.

assign values to A_-1, A₀, B_-1, and B₀, using equations 2.33 and 2.34

2.

evaluate A₁ and B₁ from the recurrence relations 2.31 and 2.32 and find A₁/B₁,

3.

evaluate A₂ and B₂ likewise, find A₂/B₂, and check for convergence, if none

4.

evaluate A₃ and B₃, as before, find A₃/B₃, and check for convergence,

5.

and so on until the result converges, i.e., until the result no longer changes, within a prescribed accuracy, as we embark on a yet another iteration.

Looking once again at the continued fraction expansion for $\Gamma (a,x)$:

$$e^{-x}x^a \left(\frac{1}{x+}\, \frac{1-a}{1+} \, \frac{1}{x+... ...}\, \cdots \, \frac{n-a}{1+} \, \frac{n}{x+}\, \cdots\right)$$

we find that

$$\begin{eqnarray*}A_{-1} &=& 1 \\ A_0 &=& b_0 = 0 \\ B_{-1} &=& 0 \\ B_0 &=& 1 \end{eqnarray*}$$

and

$$\begin{eqnarray*}a_1 &= 1 &\Rightarrow A_1 = b_1 A_0 + a_1 A_{-1} = x\cdot 0 + 1... ... n - a &\\ b_{2n} &= 1 &\\ a_{2n+1} &= n &\\ b_{2n+1} &= x \end{eqnarray*}$$

Observe that a regular sequencing begins only with a₂ and b₂, not with a₁ and b₁. In fact the first term in the continued expansion for $\Gamma (a,x)$ is the odd one. Consequently we shall begin our iterations from evaluating A₂ and B₂, and we'll use A₀, B₀, A₁, and B₁ as our startup values.

Using mathematicised notation our algorithm is going to be as follows:

$A_0 \leftarrow 0$
$B_0 \leftarrow 1$
$A_1 \leftarrow 1$
$B_1 \leftarrow x$
$R_{\textit{\scriptsize old}} \leftarrow A_1 / B_1$	assume $x \neq 0$
$n \leftarrow 1$
DO
$\qquad a_n \leftarrow n - a$	now an becomes a2
$\qquad b_n \leftarrow 1$	now bn becomes b2
$\qquad A_0 \leftarrow b_n A_1 + a_n A_0$	now A0 becomes A2
$\qquad B_0 \leftarrow b_n B_1 + a_n B_0$	now B0 becomes B2
$\qquad a_n \leftarrow n$	now an becomes a3
$\qquad b_n \leftarrow x$	now bn becomes b3
$\qquad A_1 \leftarrow b_n A_0 + a_n A_1$	now A1 becomes A3
$\qquad B_1 \leftarrow b_n B_0 + a_n B_1$	now B1 becomes B3
$\qquad R_{\textit{\scriptsize new}} \leftarrow A_1/B_1$
$\qquad \delta \leftarrow R_{\textit{\scriptsize new}} - R_{\textit{\scriptsize old}}$
$\qquad R_{\textit{\scriptsize old}} \leftarrow R_{\textit{\scriptsize new}}$
$\qquad n \leftarrow n+1$
UNTIL $\left\vert\delta/R_{\textit{\scriptsize new}}\right\vert < \epsilon$

At this stage we have carefully mapped the algorithm, without attempting any economies. But having a closer look at the DO loop as it's been outlined above shows that we can make some savings. First, we don't really need to maintain the variables a_n and b_n. Instead we can use n-a, 1, n and x directly:

DO
$\qquad A_0 \leftarrow A_1 + (n - a) A_0$	now A0 becomes A2
$\qquad B_0 \leftarrow B_1 + (n - a) B_0$	now B0 becomes B2
$\qquad A_1 \leftarrow x A_0 + n A_1$	now A1 becomes A3
$\qquad B_1 \leftarrow x B_0 + n B_1$	now B1 becomes B3
$\qquad \cdots$

This saves us 4 assignments within the loop. Assignments are not as cheap as they seem to be, because data movement is involved

Figure 2.13 shows subprogram big_gamma. You should be able to identify clearly all elements. The algorithm is exactly as specified above.

  

Figure: Subprogram big_gamma, which implements function $\Gamma (a,x)$

Having now implemented all required functions we can not only calculate the goodness of fit in our problem, which happens to be Q=0.053, but we can test the implementation of function goodness itself.

Figure 2.14 shows a short program, which calculates 1 - Q(a, x) = P(a, x), and writes the results on a data file, try.out. This file can then be displayed using Gnuplot and the graph compared against similar graphs found in the literature.

  

Figure: Program try, which evaluates P(a, x) = 1 - Q(a, x) for four different values of a: 0.5, 1.0, 3.0, and 10.0, and for $x \in [0, 15.0]$.

Figure 2.15 shows the Gnuplot command file used to convert try.out into a graph shown in Figure 2.16.

  

Figure 2.15: Gnuplot command file used to plot graph shown in Figure 2.16.

  

Figure 2.16: Function P(a, x) = 1 - Q(a, x) for a = 0.5, 1.0, 3.0, 10.0.