Difference between revisions of "Exponential Integral (a.k.a. Theis' Well Function)"
Analytophile (talk | contribs) (New page: ; 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 ...) |
Analytophile (talk | contribs) (fixed small errors (2) in c++ code) |
||
(One intermediate revision by the same user not shown) | |||
Line 18: | Line 18: | ||
double Ei(const double &x) | double Ei(const double &x) | ||
{ | { | ||
− | + | int n; | |
− | + | double Ei,term; | |
− | + | Ei=0.0; | |
− | + | if (x<15.0) | |
− | + | { | |
− | + | n=1; | |
− | + | term=-x; | |
− | + | Ei =-EULER-log(x)-term; | |
− | + | while(fabs(term)>EI_TOLERANCE) | |
− | + | { | |
− | + | term*=-x*n/(n+1)/(n+1); | |
− | + | Ei-=term; | |
− | + | n++; | |
− | + | } | |
− | + | } | |
+ | return Ei; | ||
} | } | ||
Line 60: | Line 61: | ||
===Notes=== | ===Notes=== | ||
* This algorithm implements eq. 5.1.11 from Abramowitz and Stegun's handbook of mathematical functions | * This algorithm implements eq. 5.1.11 from Abramowitz and Stegun's handbook of mathematical functions | ||
− | * The | + | * 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 | * 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 |
Latest revision as of 14:05, 3 January 2008
- 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; Ei =-EULER-log(x)-term; while(fabs(term)>EI_TOLERANCE) { term*=-x*n/(n+1)/(n+1); Ei-=term; n++; } } return Ei; }
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[edit]
- 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