Graeme Dennes has been continuing work on his Tanh-Sinh Quadrature function, and recently sent me a revised version, which can be downloaded from: Tanh-Sinh.zip
The main change is a new automatic tolerance adjustment feature which automatically increases the tolerance for badly behaved functions that will not converge to the default tolerance. The result is much better performance for a few of the sample functions. The download file also includes extensive additional notes.
Newton Excel Bach, not (just) an Excel Blog 
These routines appear to be excellent. My interest is in ground source heat exchange problems where the functions take the form of semi-infinte integrals with limits (0,∞) or (a,∞). I want to solve them on the fly using Excel spreadsheets to match field data. Running some comparisons using the exponential integral as a solution for the function f(x) = exp(-x)/x over the interval (a,∞) it appears that an upper limit of 100 appears to work okay. My questions:
(1) Are there any guidelines for setting the limits in the Tanh-Sinh integration when you are calcuating infinite or semi-infinite integrals using this excellent algorithm?
(2) The Keisan/Casio website includes a separate solution for Double Exponential Integration (a,∞) as a different algorithm from the Tanh-Sinh options. Should I be using this instead? Is it available in VBA for Excel?
Many thanks!
LikeLike
Michael – I have passed your questions onto Graham Dennes, as he is much more knowledgable on this subject than I am.
I note from the Wikipedia article (http://en.wikipedia.org/wiki/Tanh-sinh_quadrature) that Tanh-Sinh quadrature is also known as Double Exponential Integration, and maps the finite interval onto an infinite interval (-infinity, infinity). I would imagine that it would be quite straightforward to modify the method to work over a semi-infinite range, but I didn’t find a specific algorithm for doing that. In practice I don’t see any problem with simply entering a large value for the upper end of the integration range.
LikeLike
Michael and Doug:
I’m pleased to know the Tanh-Sinh integration method is of interest in your application Michael.
You’ve raised an interesting issue – the extension of the technique to semi-infinite intervals. (Also associated is the infinite (-inf,inf) interval.) As background, I personally haven’t considered such an extension to my Excel function to this time, although the 1974 paper by Takahasi and Mori develops the Double Exponential (DE) technique in principle, and discusses the transformations required for use with the intervals (a,b), (a,inf) and (-inf,inf) for mapping these to the infinite (-inf,inf) interval, then causing them to decay at a DE rate. This gives the DE technique its outstanding speed advantage. So depending upon the interval of integration required, we transform our integrand accordingly. It so happens that the Tanh-Sinh transform provides the fastest and most accurate method for integrating over the finite interval (a,b), as the 1974 paper shows. Thus the Tanh-Sinh method is but one of the DE transform methods. Other transforms are discussed in the paper for application to the other intervals, but as noted, I haven’t considered these to this time…
The Keisan/Casio site you mention (an excellent site for things mathematical) utilises the 1974 DE technique in three of its programs: the two “DE” programs, and the “Tanh Sinh” program. We simply need to choose the appropriate program depending on the interval of integration we require. By all means use these programs, as they should be mathematically correct (“perfect”), all being based on the 1974 paper. (In fact, the site uses a purpose-designed, built and patented hardware-based decimal calculating engine, so no binary, hexadecimal etc conversions or calculations are done on the way in, during, or on the way out. I know of nothing else like this available on the web.) A comparison of the results from their three programs with known, real-life functions and results would be instructive Michael.
In the meantime, if a large upper limit appears to be providing reasonable results with my Tanh-Sinh Excel program, could I suggest continuing with it, and checking it against some known functions and results. Unfortunately, I don’t have an opinion on the potential error here. Given that my program uses an automatic tolerance (and therefore error) metric, I expect the program will take care of this issue as it stands, at least by giving a result and its error estimate – if it can. Your findings here would be of great interest Michael.
I’m happy to look further into the semi-infinite and infinite interval applications for the DE method, and will report my findings in due course.
Thanks again for your interest Michael, and thank you Doug.
LikeLike
Thanks for the additional background Graeme (and sorry for getting your name wrong!).
LikeLike
Michael:
I have completed a (draft) Double Exponential integrator as a custom function for Excel, and have added it as a new worksheet in my Tanh-Sinh integrator workbook. So far it is behaving correctly. It integrates over the semi-infinite interval (a,inf), where the lower limit ‘a’ may be zero. This version cannot be used with oscillatory integrals, so it excludes integrands with eg sine and cosine functions. Like the Tanh-Sinh program, it also runs very fast.
A question. It would greatly assist me if you were able to provide me with one to several example functions (a dozen…?) of suitable integrands over the semi-infinite interval, along with their lower limits ‘a’, and their results (to as many significant digits as possible up to 15). I will use your functions as test examples to confirm proper operation, and then include them permanently in the spreadsheet as examples. If all goes well, once I add the write-up for the program’s use, source code comments, etc, which will take a few hours to finalise, I’ll seek Doug’s agreement to host the updated spreadsheet file.
Should you not have any suitable test examples, I’ll do a search through my library.
LikeLike
Gentleman,
Fantastic responses. I really appreciate your attention to and interest in this issue.
Following your advice, I have successfully examined calculation of the Exponential Integral using the Tanh_Sinh procedure as programmed by Graeme in VBA for Excel. I must say, I am simply astounded by the results!
But first, I must context this with some history. My original interest in this subject stems from practical application of the exponential integral as a basic solution to the apparently mundane question of the change in water level caused when pumping water from a well. The solution to the appropriate differential equation was originally provided in terms of heat flow by Carslaw & Jaeger in their classic tome on the Conduction of Heat.
A pioneering hydrogeologist, Charles V Theis with the USGS, figured out in 1935 that the heat flow solution for conduction in solids was a direct analogy to the groundwater flow regime created by pumping water from a well at a constant rate. His mathematical solution was called the Theis Well Function but is numerically identical to the exponential integral, and it became a cornerstone of the groundwater supply technology, especially in the analysis of pumping tests. It also got picked up by the oil industry to describe what happened in producing oil wells.
Over the next 50 or so years, tables of the function were provided at varying degrees of accuracy and hydrogeologists solved practical problems by overlaying observeed water-level data on a double-logarithmic plot of the exponential integral, in a form which was known for many years as the “Theis Type Curve“.
I should add that, over this period, the discipline of “well hydraulics” grew in leaps and bounds, developing numerous variants on the Theis Well Function to accommodate all manner of differing geologic settings hydrologic conditions and geometric situations. The challenge to using these various functions always lay in the increasingly complex differential calculus underlying the burgeoning family of solutions available.
The advent of scientific calculators and personal computers intensified the need to evaluate the well function directly, and a number of methods appeared in the literature, from polynomial approximations straight out of Abramowitz and Stegun (1965) to the simple elegance provided by the numerical inversion of the well function in Laplace space (Stehfest 1970). These proved reasonably accurate solutions but proved cumbersome — for example, in the transition from Excel 2003 to 2007, all my numerical Laplace Transform Inversion needed to be rewritten.
And so to the application of the Tanh_Sinh routine. The exponential integral is defined as the function f(x) = exp(-x)/x integrated over the interval (a,inf). I set the upper limit initially to 100 and it worked fine, except for the last value shown in this table. An upper limit of 1000 produced Excel machine accuracy (1E-15) across the entire practical range of interest. The table includes comparisons with the “perfect” calculation provided by the Double Exponential Integration routine at the Keisan.Casio website. Moreover, it is remarkably stable with an upper limit set anywhere from 1000 to 1e+9.
It would obviously be more elegant to use the Double Exponential Integration routine without having to specify some arbitrary value for infinity; I look forward to this being provided by Graeme. However, the speed, accuracy, reliability and flexibility of this method is for me simply astounding. It blows away all the other methods I have labored with for many, many years. Thank you so much for making this freely available to Excel users like me.
Regarding your final question, Graeme, I need to do some homework on the integral form of these additional functions and on rooting out the tabulated values for the functions that are of interest to me.
Many thanks.
LikeLike
Michael:
Thank you for the background and history details of your integration problem. It was indeed very informative. A real-world requirement.
Your very detailed data table comparing the results from my Tanh-Sinh function with the Keisan program was very instructive, especially regarding the accuracies achieved. Because of the high decay rate provided by the DE methods, your integrand obviously drops off very quickly, allowing a relatively small number to act as infinity, even though it is not as ideal as a true (a,inf) DE routine. Thank you for presenting your test results.
Yes, like yourself, when I achieved the initial results from the Tanh-Sinh Excel function, I could not believe what I was witnessing!! It was beyond my wildest (mathematical!) dreams. I initially thought Excel must have been drawing some data from its memory cache instead of calculating it, so I actually turned the computer off between calculation runs – just to be sure! The results remained the same.
Re the status of my current draft DE (a,inf) version. When it provides a result, the result is correct to the 15 significant digits of Excel. So far, so good. Unfortunately, for some of my test functions, it fails to provide a result (simply shows zero, but no errors), while the Keisan program does provide the correct answer for those failed test functions. I need to resolve this before I’m satisfied with it… and I can’t put a timeframe on it unfortunately. At this time, I don’t need further test functions, as all I need is some that fail – and I have a few of those!
Thanks for your interest Michael, and thank you for the positive, practical feedback on the Excel Tanh-Sinh program. I’m really pleased to know it’s been able to assist you in your work.
LikeLike
Michael:
Doug has hosted an update to my quadrature workbook. As noted there, the revised program adds two Double Exponential quadrature UDFs for the interval (a,inf). The first test integral on the worksheet DE_QUAD is your expression EXP(-x)/x.
I’d appreciate knowing if the new programs do their intended job for you.
With my thanks.
LikeLike
Pingback: Daily Download 17: Numerical Integration | Newton Excel Bach, not (just) an Excel Blog