Problem with integral on WP 34s

09122021, 09:34 PM
Post: #21




RE: Problem with integral on WP 34s
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?


09122021, 09:50 PM
Post: #22




RE: Problem with integral on WP 34s
Hi Nigel, I figured out how to check the flag setting with the FS? option. My Flag D is set successfully. Thanks so much!


09122021, 11:15 PM
(This post was last modified: 09132021 09:39 AM by Albert Chan.)
Post: #23




RE: Problem with integral on WP 34s
(09122021 08:41 PM)Nigel (UK) Wrote: The comments from [double exponential] integration code state that Even if integrand is analytic, it may be wise to breakup 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, breakup 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.0e24')), 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') 

09132021, 10:20 AM
Post: #24




RE: Problem with integral on WP 34s
You can find a version of the WP34S for DM42 that includes Romberg integration here:
https://gitlab.com/njdowrick/wp34sford...ee/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/389...ndows/bin/ This has the Romberg routine as its standard integration technique, and it will be far faster than either physical calculator. Nigel (UK) 

09152021, 02:14 AM
Post: #25




RE: Problem with integral on WP 34s
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! 

09152021, 02:21 AM
Post: #26




RE: Problem with integral on WP 34s
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!)?


09152021, 07:41 AM
Post: #27




RE: Problem with integral on WP 34s
Hi, lrdheat
For tanhsinh, a to b, c = (a+b)/2, sample points f(c), then f(cd*r), f(c+d*r) For sinhsinh, ∞ to ∞, c = 0, sample points f(c), then f(cd*r), f(c+d*r) For expsinh, a to ∞, c = a, sample points f(c+d), then f(c+d/r), f(c+d*r) This may be oversimplified. see double_exponential.py https://www.hpmuseum.org/forum/thread16549.html https://www.hpmuseum.org/forum/thread16...#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.87E217142, 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.7892954451426135e8'), 93, 4, 0, [ mpf('2.458125108263046'), mpf('1.2794193804039591'), mpf('1.0055096455855115'), mpf('1.0000000178929545'), mpf('1.0') ]) 

09152021, 11:26 AM
Post: #28




RE: Problem with integral on WP 34s
(09152021 02:21 AM)lrdheat Wrote: How, then, does the WP 34S on DM 42 recognize the important part of the curve when Quote:For expsinh, 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, expsinh 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.5934163150981393e009 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.448208716454109e013 61 

09162021, 12:27 PM
Post: #29




RE: Problem with integral on WP 34s
(09152021 02:14 AM)lrdheat Wrote: Hi again, Nigel, For the doubleexponential routine, if you interrupt by pressing < (backspace) the most recent estimate is left in the xregister. 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. (09152021 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 doubleexponential integration method is, and how well it is implemented on the WP34S. 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 nonanalytic) somewhere in the range of integration. I'm happy to spot these points myself and to work around them. Nigel (UK) 

09162021, 04:54 PM
Post: #30




RE: Problem with integral on WP 34s
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 WP34S DM 42!


09172021, 10:02 AM
Post: #31




RE: Problem with integral on WP 34s
(09162021 04:54 PM)lrdheat Wrote: [...] a nicely implemented integration technique! (09162021 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 WP34S I am the one that wrote the wp34s doubleexponential 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. César  Information must flow. 

09182021, 11:38 AM
Post: #32




RE: Problem with integral on WP 34s
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) 

09182021, 03:53 PM
Post: #33




RE: Problem with integral on WP 34s
An example of how nicely this implementation works on the DM 42 version of WP34S is the quick convergence to pi in the integral of 1/(sqrt x * (1+x)) from 0 to +inf.
On my CASIO fx9750iii, 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+. 

09182021, 04:38 PM
Post: #34




RE: Problem with integral on WP 34s
Furthermore, the CASIO times out without calculating the integral if the upper limit is >1842. At 1842, the integral is 3.095


09182021, 04:45 PM
(This post was last modified: 09182021 05:35 PM by Albert Chan.)
Post: #35




RE: Problem with integral on WP 34s
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_infiniteintegrals.html 

09182021, 05:11 PM
Post: #36




RE: Problem with integral on WP 34s
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 WP34S handles the 0 to inf domain in quick, excellent fashion! 

09182021, 05:16 PM
Post: #37




RE: Problem with integral on WP 34s
…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!


09182021, 06:26 PM
Post: #38




RE: Problem with integral on WP 34s
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 7E33).


« Next Oldest  Next Newest »

User(s) browsing this thread: