[wp34s] New integration program
|
03-08-2017, 05:46 PM
(This post was last modified: 03-13-2017 05:59 AM by emece67.)
Post: #1
|
|||
|
|||
[wp34s] New integration program
Hi all,
In another thread I presented a description of an algorithm (not mine) to approximate definite integrals different to the Romberg method classically used by hp machines: the tanh-sinh method. After some time I was able to write a program for the wp34s machine that uses this method to compute definite integrals. Find it below. In my tests this method is much faster than the Romberg one & seems to manage appropriately a wider class of integrands (infinite integrands at one or both interval edges, infinite derivatives at one or both interval ends). I really think that this method may replace the Romberg one even in small handheld calculators, but it is the first implementation of such method I know in a calculator & much more tests are needed to check its behaviour, so beta testers are welcome. Regards. Code:
NOTE: updated to v1.0l (some tweaks in order to speed-up execution, get more precision in some difficult cases and a more responsive user interrupt) Note 2: updated to work in DBLON mode. |
|||
03-10-2017, 01:35 PM
(This post was last modified: 03-12-2017 10:43 PM by emece67.)
Post: #2
|
|||
|
|||
RE: [wp34s] New integration program
Find attached my findings when comparing this new program against the built in integration program.
Tanh-Sinh_comparison.pdf (Size: 21.23 KB / Downloads: 72) Finally I used a battery of 40 test integrals of various kinds. On this table:
When counting the correct digits I was only concerned about the 12 visible digits (if the expected result is 0, the table shows the absolute value of the exponent of the returned result). Time was measured with a handheld chrono, so don't expect high accuracy on such. Times marked as ">xxx" means that I stopped measuring the time when xxx seconds get reached. "est. time" means that I estimated the time required for the built-in integrator by measuring a few iterations on the real machine and using the simulator to get the number of total iterations needed. The machine was always in SLOW and DBLOFF modes. Regards. NOTE: updated to show results of last v1.0l version. Some typos fixed. |
|||
03-11-2017, 12:22 PM
Post: #3
|
|||
|
|||
RE: [wp34s] New integration program
Any chance of adding the WP43 version 2.2 firmware to the results? The integration there was a non-adaptive Gauss-Kronrod 10/21 point quadrature -- it is exact for polynomials up to degree 19, although subject to numeric stability, but not so great for other functions.
Pauli |
|||
03-11-2017, 12:45 PM
Post: #4
|
|||
|
|||
RE: [wp34s] New integration program
(03-11-2017 12:22 PM)Paul Dale Wrote: Any chance of adding the WP43 version 2.2 firmware to the results? So WP43 firmware is already at version 2.2?!? Great news! Just kidding, SCNR. Greetings, Massimo -+×÷ ↔ left is right and right is wrong |
|||
03-11-2017, 03:10 PM
Post: #5
|
|||
|
|||
RE: [wp34s] New integration program
V2.2 used the 10/21 Gauss-Kronrod method. Current v3.3 uses the Romberg Method. The table I posted compares my implementation of the Tanh-Sinh method against the Romberg method in v3.3.
|
|||
03-12-2017, 11:07 PM
(This post was last modified: 03-13-2017 06:03 AM by emece67.)
Post: #6
|
|||
|
|||
RE: [wp34s] New integration program
I've just updated the program listing and the comparison table cause I've added some tweaks to the program.
Overall I think that it, definitely, wins against the built-in, Romberg-based integration program in the v3.3 wp34. The main source of trouble for this algorithm seems to be such functions with infinite derivatives inside the integration interval (even if the function does not go to infinity). There are other cases where the built-in integrator is better, but even in such cases this tanh-sinh method gives acceptable results (F22) or, simply, both methods are failing miserably (F30, F32). In all these problem cases, switching to DBLON mode does not help. On the other side, the huge speed-up in almost any situation of this tanh-sinh method and its better results when the integrand goes to infinity at the integration interval edges makes it, IMHO, a safe replacement for the built-in integrator in the wp34 calculator and, perhaps, in other machines with Romberg based integrators. Even more, with minor changes, this method can cope with integrals where one of the limits is infinite (exp-sinh method) or both are (sinh-sinh method). Regards. |
|||
03-17-2017, 06:16 PM
Post: #7
|
|||
|
|||
RE: [wp34s] New integration program
I've continued looking for ways to improve this program. Although I've not yet implemented it for the wp34s, I do now have an emulation, in Python, of a new version of the program that also manages not only finite integration intervals, but both infinite and semi-infinite ones.
I'll post the wp34s version in a few days. Regards. |
|||
03-24-2017, 05:12 PM
(This post was last modified: 03-27-2017 05:52 PM by emece67.)
Post: #8
|
|||
|
|||
RE: [wp34s] New integration program
I've just posted an implementation of the double-exponential method (that is: tanh-sinh, exp-sinh or sinh-sinh, depending on the integration interval ends being finite or not) in the General Software library.
This new program has grown in size (about 50 steps for a total of ~250) but now it copes with infinite or semi-infinite integration intervals. It can now manage things like: \[\int_{-1}^\infty{\ln(x+1)\over x^2+1}dx={3\pi\ln(2)\over8}\approx0.8165947838638507989377583368391052\ldots\] that has a discontinuity at one interval end (-1), being the other end infinite, getting as result:
Although there are, of course, integrands that may drive this method in trouble, I'm amazed about its performance and capabilities. Regards. |
|||
08-11-2019, 05:33 AM
Post: #9
|
|||
|
|||
RE: [wp34s] New integration program
(03-17-2017 06:16 PM)emece67 Wrote: I've continued looking for ways to improve this program. Although I've not yet implemented it for the wp34s, I do now have an emulation, in Python, of a new version of the program that also manages not only finite integration intervals, but both infinite and semi-infinite ones. Yes, it's been a while since the above message was placed! The post discusses a new Python version which manages semi-infinite and infinite intervals as well as finite intervals. Is the Python code listing available? Thanks nova |
|||
04-03-2021, 01:40 PM
Post: #10
|
|||
|
|||
RE: [wp34s] New integration program
(08-11-2019 05:33 AM)nova Wrote: Is the Python code listing available? Hi, As I have stated in another thread, my original Python code for the double-exponential quadrature is now on GitHub. It has suffered some refactoring from its original shape (that I used to derive the wp34s code). Regards. |
|||
04-04-2021, 03:33 PM
(This post was last modified: 04-10-2021 10:45 PM by Albert Chan.)
Post: #11
|
|||
|
|||
RE: [wp34s] New integration program
(03-24-2017 05:12 PM)emece67 Wrote: \[\int_{-1}^\infty{\ln(x+1)\over x^2+1}dx={3\pi\ln(2)\over8}\approx0.8165947838638507989377583368391052\ldots\] With this integral transformed, I noted some weekness of tanh-sinh. (Or, it might be just a fluke with this example ...) Rewrite above integral, then let x=e^y: \(\displaystyle \int _0 ^∞ {\ln(x) \over 1+(x-1)^2} dx = \int _{-∞} ^∞ {y\;e^y \over 1+(e^y-1)^2} dy\) >>> from mpmath import * >>> import double_exponential # emece67 source from here >>> quadde = lambda f,(a,b): double_exponential.double_exponential(f,a,b) >>> f = lambda y: y*exp(y) / (1 + expm1(y)**2) >>> quadde(f, [-inf,inf]) # method sinh-sinh (mpf('0.81659478386382767'), mpf('7.7731406533665393e-8'), 73, 4, 2) To use tanh-sinh, just convert back intervals: [-∞,+∞] -> [0,1] >>> def f01(t): y=t*t-t; return f((2*t-1)/y) * (2*y+1)/(y*y) ... >>> quadde(f01, [0,1]) # method tanh-sinh (mpf('0.81659478390520546'), mpf('1.0424940875997102e-5'), 77, 5, 0) It seems sinh-sinh does better than tanh-sinh, for interval [-∞,+∞] --- >>> quadde(f, [-40,40]) (0, mpf('0.83732598514658618'), 87, 5, 0) A bit unexpected, getting area of 0.0, with huge error estimate. plot(f, [-40,40]), showed it has a shape of a heartbeat (inside ±5) We should instead move dominant part to the edge. >>> fold = lambda x: f(x) + f(-x) >>> quadde(fold, [0,40]) (mpf('0.81659478386385076'), mpf('2.0486468077507425e-10'), 161, 5, 0) Result is similar to exp-sinh, for [0,+∞] >>> quadde(fold, [0,inf]) (mpf('0.81659478386385087'), mpf('4.1397899552819695e-9'), 131, 4, 1) Edit: I misinterpreted last argument, should be 0=tanh-sinh, 1=exp-sinh, 2=sinh-sinh |
|||
04-04-2021, 05:09 PM
(This post was last modified: 04-04-2021 05:16 PM by emece67.)
Post: #12
|
|||
|
|||
RE: [wp34s] New integration program
(04-04-2021 03:33 PM)Albert Chan Wrote: To use tanh-sinh, just convert back intervals: [-∞,+∞] -> [0,1] I suppose it depends on the particular change of variable you use to get from [-∞, +∞] to [0, 1]. Maybe you can find cases when the tanh-sinh method behaves better in [0, 1] than transforming the problem to the exp-sinh case. (04-04-2021 03:33 PM)Albert Chan Wrote: >>> quadde(f, [-40,40]) This implementation of the double-exponential algorithm is aimed at a handheld calculator, so it uses simple heuristics in some places. Here you are fooling one of such heuristics (related to the fact that the function vanishes rapidly away from 0). The first iterations of the algorithm "doesn't see" the interesting part of the function, and when later iterations find something, the algorithm thinks it is round-off noise. In some way this function is "bell-shaped", and this algorithm concentrates the sampling points at the extremes of the integration interval, so caution is advised for cases such this one. Integrating exp(-x^2) in [-10, 10] or [-inf, +inf] gives good results, but in [-20, 20] triggers the heuristic. You have found a way to circumvent this, though. Regards. |
|||
« Next Oldest | Next Newest »
|
User(s) browsing this thread: