+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: Not HP Calculators (/forum-7.html)
+--- Forum: Not quite HP Calculators - but related (/forum-8.html)
+--- Thread: Problem with integral on WP 34s (/thread-17448.html)
Pages: 1 2
RE: Problem with integral on WP 34s - lrdheat - 09-12-2021 09:34 PM
Thanks! Just to make sure that I am successful in setting flag D on the DM 42 version of WP 34S, how do I enter “D” on the SF prompt?
RE: Problem with integral on WP 34s - lrdheat - 09-12-2021 09:50 PM
Hi Nigel, I figured out how to check the flag setting with the FS? option. My Flag D is set successfully. Thanks so much!
RE: Problem with integral on WP 34s - Albert Chan - 09-12-2021 11:15 PM
(09-12-2021 08:41 PM)Nigel (UK) Wrote: The comments from [double exponential] integration code state that
Quote:the double exponential method relies on the function to be
integrated being analytic over the integration interval, except,
perhaps, at the interval ends
Even if integrand is analytic, it may be wise to break-up integral.
Example, ∫(z^(1/3), z = -2 .. 3)
>>> from mpmath import *
>>> def f(x): f.n+=1; return cbrt(x) # complex cube root, principle branch
...
>>> f.n=0; quad(f, (-2,3), error=1), f.n
((mpc(real='4.1920063238379486', imag='1.6347100863481141'), mpf('0.01')), 427)
Now, break-up integral to 2 pieces.
>>> f.n=0; quad(f, (-2,0,3), error=1), f.n
((mpc(real='4.1900023206128241', imag='1.6366854539575821'), mpf('1.0e-24')), 106)
>>> F = lambda z: mpc(z)**(4/3)/(4/3)
>>> F(3) - F(-2) # integral true value
mpc(real='4.1900023206128232', imag='1.6366854539575817')
>>> F(3) + F(2)*exp(pi/3*1j) # another way
mpc(real='4.1900023206128232', imag='1.6366854539575821')
RE: Problem with integral on WP 34s - Nigel (UK) - 09-13-2021 10:20 AM
You can find a version of the WP34S for DM42 that includes Romberg integration here:
https://gitlab.com/njdowrick/wp34s-for-dm42/-/tree/romberg
I've made it a separate branch rather than including it in the Master branch to keep the main code relatively uncluttered. The command is called \({\rm ROM}\int\) and it is in the Integrate menu (shift 8).
I'm not sure how useful this will be - although the DM42 is faster than the WP34S when I've compared them, this integral still seems to be taking a long time to converge even with the Romberg routine. If you don't want to wait, here's a link to the folder of the v3892 WP34S Sourceforge repository where the Windows emulator lives:
https://sourceforge.net/p/wp34s/code/3892/tree/trunk/windows/bin/
This has the Romberg routine as its standard integration technique, and it will be far faster than either physical calculator.
Nigel (UK)
RE: Problem with integral on WP 34s - lrdheat - 09-15-2021 02:14 AM
Hi again, Nigel,
If an integral is taking a fair amount of time to run, and I am satisfied that the estimate on screen is a good one, is there a way to retrieve the latest estimate? When I stop with r/s or exit, the estimate must reside somewhere!
RE: Problem with integral on WP 34s - lrdheat - 09-15-2021 02:21 AM
Unrelated, but interesting…integral of x*e^-x from 0 to infinity should be 1. When integrating from, say 0 to 1,000, it converges on 1. If I integrate from 0 to 1,000,000, it converges on 0, misses the important part of the curve close to 0. How, then, does the WP 34S on DM 42 recognize the important part of the curve when integrating from 0 to infinity (it converges on 1!)?
RE: Problem with integral on WP 34s - Albert Chan - 09-15-2021 07:41 AM
Hi, lrdheat
For tanh-sinh, a to b, c = (a+b)/2, sample points f(c), then f(c-d*r), f(c+d*r)
For sinh-sinh, -∞ to ∞, c = 0, sample points f(c), then f(c-d*r), f(c+d*r)
For exp-sinh, a to ∞, c = a, sample points f(c+d), then f(c+d/r), f(c+d*r)
This may be over-simplified. see double_exponential.py
https://www.hpmuseum.org/forum/thread-16549.html
https://www.hpmuseum.org/forum/thread-16635-post-148080.html#pid148080
f(x) = x*exp(-x), important part is from 0 to 10, peak at 1.0
After peak, it decay fast
f(0) = 0
f(1) ≈ 0.368 (peak)
f(10) ≈ 0.000454
f(20) ≈ 0.0000000412
f(30) ≈ 0.00000000000281
Convergence is based from previous estimate, with 1st sample point f(c)
If we integrate from 0 .. 1e6, f(c=5e5) ≈ 2.87E-217142, underflow to 0.0
Next estimate also 0.0 ... it "converged"
Or course, if system designed not to underflow, it still work (but not great)
With more reasonable limits, we get better result.
>>> from double_exponential import *
>>> f = lambda x: x*exp(-x)
>>> double_exponential(f, 0, 1e6)
(mpf('0.99999999958767749'), mpf('0.00028631970183223832'), 197, 5, 0,
[
mpf('0.019275511359964845'),
mpf('0.009934112577158722'),
mpf('0.40518132830865922'),
mpf('1.0876426215773149'),
mpf('0.99971367988584525'),
mpf('0.99999999958767749')
])
>>> double_exponential(f, 0, 100)
(mpf('1.0'), mpf('1.7892954451426135e-8'), 93, 4, 0,
[
mpf('2.458125108263046'),
mpf('1.2794193804039591'),
mpf('1.0055096455855115'),
mpf('1.0000000178929545'),
mpf('1.0')
])
RE: Problem with integral on WP 34s - Albert Chan - 09-15-2021 11:26 AM
(09-15-2021 02:21 AM)lrdheat Wrote: How, then, does the WP 34S on DM 42 recognize the important part of the curve when
integrating from 0 to infinity (it converges on 1!)?
Quote:For exp-sinh, a to ∞, c = a, sample points f(c+d), then f(c+d/r), f(c+d*r)
Sorry, I missed the question.
r → ∞: (c+d/r) is approaching c, (c+d*r) to infinity. (for a to ∞, d>0)
r → 1 : (c+d/r) and (c+d*r) are both approaching (c+d)
So, exp-sinh is really summing of 2 integrals.
\( \displaystyle \int_c^∞ = \int_c^{c+d} + \int_{c+d}^∞ \)
Test using code from here (I had updated quad to return function call counts too)
lua> Q = require'quad'
lua> f = function(x) return x*exp(-x) end
For x = 0 .. inf, f(x) peaked at x=1, so default d=1 is not too bad.
lua> Q.quad(f, 0, huge) -- d = 1 (default)
0.9999999999749527 2.5934163150981393e-009 107
But, from plots, we know 0 .. 10 is where the action is.
lua> Q.quad(f, 0, huge, 10) -- d = 10
0.999999999999988 7.448208716454109e-013 61
RE: Problem with integral on WP 34s - Nigel (UK) - 09-16-2021 12:27 PM
(09-15-2021 02:14 AM)lrdheat Wrote: Hi again, Nigel,
If an integral is taking a fair amount of time to run, and I am satisfied that the estimate on screen is a good one, is there a way to retrieve the latest estimate? When I stop with r/s or exit, the estimate must reside somewhere!
For the double-exponential routine, if you interrupt by pressing <- (backspace) the most recent estimate is left in the x-register. There doesn't seem to be anything similar in the Romberg routine. These routines use local registers and flags, so their internal workings are only accessible if the routine allows them to be.
(09-15-2021 02:21 AM)lrdheat Wrote: Unrelated, but interesting…integral of x*e^-x from 0 to infinity should be 1. When integrating from, say 0 to 1,000, it converges on 1. If I integrate from 0 to 1,000,000, it converges on 0, misses the important part of the curve close to 0. How, then, does the WP 34S on DM 42 recognize the important part of the curve when integrating from 0 to infinity (it converges on 1!)?
When both limits are infinite, the integration routine does a change of variable from \(x\) to \(t\), where
\[x=\sinh\left((\pi/2)\sinh(t)\right)\]
(or something similar). When \(t=1\), \(x=3.1\); when \(t=5\), \(x=2.1\times10^{50}\)! So it's only necessary to integrate over a small range of \(t\) to reach essentially infinite values of \(x\). This means that the important region for \(x\) will automatically be found, so long as it's concentrated somewhere near \(x=0\) as it is in this case.
The more I've read about this subject, the more I am impressed by how good the double-exponential integration method is, and how well it is implemented on the WP-34S. I won't be using Romberg myself - the cases for which the double exponential method fails always seem to be those where the function being integrated or its derivatives are discontinuous (i.e., the function is non-analytic) somewhere in the range of integration. I'm happy to spot these points myself and to work around them.
Nigel (UK)
RE: Problem with integral on WP 34s - lrdheat - 09-16-2021 04:54 PM
Thanks Albert and Nigel for the excellent explanations of a nicely implemented integration technique! I agree that it is more efficient and just about as easy to break the integral up into multiple integrals when encountering discontinuities! Glad that this was included in the WP-34S DM 42!
RE: Problem with integral on WP 34s - emece67 - 09-17-2021 10:02 AM
(09-16-2021 04:54 PM)lrdheat Wrote: [...] a nicely implemented integration technique!
(09-16-2021 12:27 PM)Nigel (UK) Wrote: The more I've read about this subject, the more I am impressed by [...] how well it is implemented on the WP-34S
I am the one that wrote the wp34s double-exponential routine, so I'm really pleased by your comments and happy that you find it useful.
A few months ago there was a thread here where people (robve and Albert Chan) discussed ways to improve the algorithm. I'm no longer in position to continue developing the routine (in fact I removed the wp34s toolchain from mi computer months ago, simply lost interest...) but from the comments on such thread, there is room for improvement, so I encourage any interested on it to read the previous thread.
Regards.
RE: Problem with integral on WP 34s - Nigel (UK) - 09-18-2021 11:38 AM
An excellent thread which I completely missed at the time. Thank you for the link. I may have a go at implementing some of the (modest) improvements that have been suggested, once I understand it all!
Nigel (UK)
RE: Problem with integral on WP 34s - lrdheat - 09-18-2021 03:53 PM
An example of how nicely this implementation works on the DM 42 version of WP-34S is the quick convergence to pi in the integral of 1/(sqrt x * (1+x)) from 0 to +inf.
On my CASIO fx-9750iii, even with a domain as wide as from 0 to 1000, after ~6 seconds, the integral is still well shy of pi, coming in with 3.078+.
RE: Problem with integral on WP 34s - lrdheat - 09-18-2021 04:38 PM
Furthermore, the CASIO times out without calculating the integral if the upper limit is >1842. At 1842, the integral is 3.095
RE: Problem with integral on WP 34s - Albert Chan - 09-18-2021 04:45 PM
CASIO is correct !
Let y = x*x, dy = 2x dx, to remove √(x)
∫(1/(√(x)*(1+x)), x) = ∫(2/(1+y*y), y) = 2*atan(y)
I(x = 0 .. 1000) = 2*atan(√1000) ≈ 2*(pi/2 - 1/√1000) ≈ 3.078
I(x = 0 .. 1842) ≈ 2*atan(√1842) ≈ 2*(pi/2 - 1/√1842) ≈ 3.095
To numerically calculate the integral (either x or y version), noted that \(\int_0^a = \int_{1/a}^∞ \quad\), a ≥ 0
For calculators without 'inf' feature, we can transform it back to 0 .. 1
http://fmnt.info/blog/20180818_infinite-integrals.html
RE: Problem with integral on WP 34s - lrdheat - 09-18-2021 05:11 PM
CASIO is correct for that domain. The CASIO times out for an upper limit if >1842 with a 3.095 calculation. If I then do the integral over a domain of 1843 to 200000,
the 2 integrals added together are still shy, at ~3.135+. The WP-34S handles the 0 to inf domain in quick, excellent fashion!
RE: Problem with integral on WP 34s - lrdheat - 09-18-2021 05:16 PM
…even using a domain of 1843 to 10000000, The 2 integrals added do not get me to a correct approximation of pi to the 4 th digit after the decimal point!
RE: Problem with integral on WP 34s - lrdheat - 09-18-2021 06:26 PM
Meant to break it into 2 integrals, from 0 to 1842 + from 1842 to 10000000. This yields ~3.14096. The WP 34S on DM 42 from 0 to inf shows 3.141592653+ (differed from pi by 7E-33).