Exponential Integral (a.k.a. Theis' Well Function)

From aemwiki
Revision as of 08:03, 21 December 2007 by Analytophile (talk | contribs) (Notes)

Jump to: navigation, search
Language
C++, Visual Basic
Purpose
Calculate the Exponential integral function, <math>\mathbf{E_1}(x)</math>, also known as the Theis Well function <math>W(u)</math>
Authors
Theo Olsthoorn (VB), James R. Craig (c++)

C++ Version:

#include <math.h>
const double EI_TOLERANCE=1e-9;
const double EULER=0.57721566490153;
/***********************************************************************
	 Exponential Integral / Well Function
************************************************************************
Exponential Integral, E_1(x)=Real(-Ei(-x))
Author: James R. Craig
translated from the VB version published in T. Olsthoorn, "Do a bit more 
with convolution", Groundwater 46(1), 2008
-----------------------------------------------------------------------*/
double Ei(const double &x)
{ 
  int    n;
  double Ei,term;
  Ei=0.0;
  if (x<15.0) 
  {
    n=1;
    term=-x;
    W =-EULER-log(x)-term;
    while(fabs(term)>EI_TOLERANCE)
    {
      term=-u*n/(n+1)*n/(n+1)* term;
      Ei-=term;
      n++;
    }
  }
} 

Visual Basic Version:

'***********************************************************************
'	 Theis Well Function
'************************************************************************
'Exponential Integral
'from T. Olsthoorn, "Do a bit more with convolution", Groundwater 46(1), 2008
'-----------------------------------------------------------------------*/
Function W (u As Double) As Double
Dim n AS Integer, term As Double
W = 0
If u < 15 Then
 n = 1
 term = -u
 W =-0.5772156649 - Log(u) - term
 Do While Abs(term) > 0.000000001
  term = -u * n/(n + 1) ^ 2 * term
  W = W- term
  n = n + 1
  Loop
 End If
End Function

Notes

  • This algorithm implements eq. 5.1.11 from Abramowitz and Stegun's handbook of mathematical functions
  • The approximation used yields an error bound of <math>10^{-9}</math>
  • These implementations do not address some important issues (i.e., x<0 shouldn't be allowed). A more robust implementation is available from Press et. al, Numerical Recipes