HP71B Integral Questions

02022020, 03:31 PM
Post: #1




HP71B Integral Questions
Tried on Emu71b, energy integral from thread Sharp EL506 vs Sharo EL516T
From the user's manual, about the algoritm, it uses Romberg's method, after a nonlinear transform. Code: 10 A=6371000 @ B=9.46073047258E+15 @ C=3.98589196E+17 @ A1=A1 >RUN P=?1E5 IA( 127 ) = 62563050744 625265.300862 .11 SEC IB( 255 ) = 62563050665 625640.053047 .3 SEC >RUN P=?1E6 IA( 127 ) = 62563050647 62526.5300862 .11 SEC IB( 255 ) = 62563050665 62564.0053047 .3 SEC >RUN P=?1E7 IA( 255 ) = 62563050656 6255.39238586 .22 SEC IB( 255 ) = 62563050665 6256.40053047 .31 SEC Why does IA(127) have different area, for P=1E5 vs P=1E6 ? Both are using exactly the same 127 sample points. Does IBOUND (second column) has something to do with different IA(127) results ? How is IBOUND intepreted ? At what IBOUND numbers should an integral be splitted ? 

02032020, 01:58 PM
(This post was last modified: 02032020 02:29 PM by Albert Chan.)
Post: #2




RE: HP71B Integral Questions
(02022020 03:31 PM)Albert Chan Wrote: Why does IA(127) have different area, for P=1E5 vs P=1E6 ? Answering myself, first column = trapezoidal area of utransformed FNA(X) Code: Intervals Richardson Extrapolations Note: utransformed integrand returned 0 for integrand limit of 1 and 1, thus 1 trapezoid have zero area. Extending rows of 128 intervals, extrapolate to its limits, we have: 62522713769 Δ 62563113307 40399538 62563049732 63575 62563050740 1007 62563050648 91 62563050583 65 <\ 62563050564 19 < extrapolations does not improve estimate, actually made it worse 62563050559 5 </ The numbers are done in Excel, which use IEEE754 binary64. Except for slight rounding errors, the bold numbers matched IA(127) for P=1E5 and 1E6. 

02032020, 03:00 PM
Post: #3




RE: HP71B Integral Questions
(02032020 01:58 PM)Albert Chan Wrote: Answering myself, first column = trapezoidal area of utransformed FNA(X) I'm interested in documenting the Math ROM algorithms, can you please post the Excel formulae you used to reconstruct the INTEGRAL results? From the HP Journal article, July 1984 p31 on the Math ROM "A Closer Look at INTEGRAL": "INTEGRAL proceeds by computing a sequence of weighted sums [...]. These sums are accumulated in an extendedprecision, eightlevel Romberg scheme.[...] These iterations produce a sequence of approximate integrals I(k,j), where k=0,1,... and j=0,1,...,min(k,7). Here, I(k,0) is the direct weighted sum of 2^k samples integration points, and I(k,j) is the Romberg extrapolation given by I(k,j) = (4^j.I(k,j1)  I(k1,j1)) / (4^j1)." Is the 8level Romberg scheme close to the Richardson Extrapolations you used? About IBOUND: IBOUND is an estimation of the absolute error of the computed integral. Some details on how it is computed are also mentioned in the article cited above. IBOUND is usually close to the userprovided target precision P times the integral estimation. JF 

02032020, 04:15 PM
(This post was last modified: 04142021 12:52 AM by Albert Chan.)
Post: #4




RE: HP71B Integral Questions
(02032020 03:00 PM)JF Garnier Wrote: Is the 8level Romberg scheme close to the Richardson Extrapolations you used? Romberg's method ≡ Richardson Extrapolations applied to Trapezoids area. I did the utransformed trapezoids area in Lua, then copy/paste to Excel Code: function u(f, a, b)  utransformed f lua> a, b, c = 6.371e6, 9.4607304725808e15, 3.98589196e17 lua> exp, log = math.exp, math.log lua> function fna(x) return c/exp(x) end lua> t = trapezoid(u(fna, log(a), log(b)), 1, 1) lua> for i=1,7 do print(2^i, t()) end 2 25715027.086919405 4 13721555535.80971 8 50890703123.90066 16 59914609212.0083 32 61913859849.80031 64 62401515156.47649 128 62522713769.32443 Only 2 formulas used in Excel, one for the weight, second for the extrapolation. Others are just copies. With this setup, Simpson's rule have weight of 3, Boole's rule have weight of 15, ... Code: A B C 

02032020, 11:15 PM
(This post was last modified: 02032020 11:40 PM by Albert Chan.)
Post: #5




RE: HP71B Integral Questions
(02022020 03:31 PM)Albert Chan Wrote: From the user's manual, about the algoritm, it uses Romberg's method, after a nonlinear transform. There was an error in the manual. x = ½(3uu³) substitution *required* integrand limit 1 to 1, not a to b \(\large \int _{1} ^1 f(x) dx = {3 \over 2}\int _{1} ^ 1 (1u^2) f \left( { u (3u^2) \over 2} \right) du \) For a general case, with integrand limit a to b, we have: \(\large \int _a ^b f(x) dx = {3(ba) \over 4}\int _{1} ^ 1 (1u^2) f \left( { 2(b+a) + (ba) u (3u^2) \over 4} \right) du \) From the manual: Quote:INTEGRAL uses extended precision. Internally, sums are accumulated in 15digit numbers ... Trivia: For 2 intervals, h = 1 → 2^16 intervals, h = 1/2^15 = 0.00003 05175 78125 (exact) With 15 digits precision, and u ≤ 1, this implied INTEGRAL internally calculated u's are all exact. 

02052020, 08:43 AM
Post: #6




RE: HP71B Integral Questions
(02032020 11:15 PM)Albert Chan Wrote:(02022020 03:31 PM)Albert Chan Wrote: From the user's manual, about the algoritm, it uses Romberg's method, after a nonlinear transform. Thanks, Albert. Actually, this description originates from a HP15C manual (text originally written by W. Kahan) ... and was simplified in the HP71 Math ROM manual. In the HP15C Advanced Functions Handbook p.41, the integration limits are 1 and 1, and the exact sentence is: "Their spacing can be demonstrated by substituting, say" (emphasis added) to indicate it is an example for illustration only, not an exact detailed description, The important word "say" was removed in the HP75 and HP71 manuals... Quote:Trivia:Interesting. The older HP34C and HP15C (and HP41C Advantage ROM) used 13 digits internally, maybe were they more limited in the number of intervals. JF 

02052020, 05:09 PM
(This post was last modified: 02092020 11:31 PM by Albert Chan.)
Post: #7




RE: HP71B Integral Questions
I noticed utransformation can be applied more than once, accelerate some integrals.
Example: \( \int _{1} ^1 {dx \over \sqrt{1x^2}} = 2\;\sin^{1} 1 = \pi\) Below, FNB(U) is utransformed FNA(X), which is more accurate, and run much faster. Code: 10 INPUT "P=? ";P >RUN P=? 1E4 IA( 8191 ) = 3.14133371246 3.14116984016E4 6.34 SEC IB( 31 ) = 3.14159247755 3.14127673304E4 .03 SEC >RUN P=? 1E5 IA( 65535 ) = 3.14156045534 3.1415397954E5 50.97 SEC IB( 31 ) = 3.14159263336 3.14127673304E4 .03 SEC Lets call integrand of IB (= utransformed FNB) as FNC. If we plot FNB and FNC, u = 1 to 1, FNB is much flatter, min=FNB(0)=3/2, max=FNB(±1)=√3 Intuition may suggest IA is easier to integrate, but it will be dead wrong. The code never evaluate the endpoints integrand limit, and skip over them. That's why we only have 2^N1 samples points, instead of 2^N+1. → Trapezoid number, T_{1} = 0 (*always*) FNB, which T_{1} should really be 2√3, skewed badly all the romberg numbers. To confirm, I tried doing by hand, with T_{1}=2√3, IA(31+2) = 3.14159265428 FNC does not suffer from this, since limit(FNC, U=1) = limit(FNC, U=+1) = 0 

02052020, 06:20 PM
Post: #8




RE: HP71B Integral Questions
I ran the same tests on the HP75 Math ROM, here is the program adapted from Albert's code for the 75C:
10 A=6371000 @ B=9.46073047258E15 @ C=3.98589196E17 @ A1=A1 20 DEF FNA(X) 22 N=N+1 @ FNA=C/EXP(X) 25 END DEF 30 DEF FNB(X) 32 N=N+1 @ X=EXP(X) @ FNB=C*X/(X+A1)^2 35 END DEF 40 INPUT "P=?";P 50 T=TIME @ N=0 @ I=INTEGRAL(LOG(A),LOG(B),P,FNA(X)) 60 DISP "IA(";N;") =";ROUND(I,0),IBOUND,TIMET;"SEC" 70 ! T=TIME @ N=0 @ I=INTEGRAL(0,LOG(BA1),P,FNB(X)) 80 ! DISP "IB(";N;") =";ROUND(I,0),IBOUND,TIMET;"SEC" The results are a bit different from the HP71 outputs: P=?2e2 IA( 63 ) = 62563132323 422758135.629 .11 SEC P=?1e3 IA( 127 ) = 62563050559 21501195.6588 .22 SEC P=?1e4 IA( 255 ) = 62563050656 2072473.90909 .495 SEC P=?1e5 IA( 255 ) = 62563050656 207247.390909 .494 SEC P=?1e6 IA( 511 ) = 62563050656 20718.6780627 .989 SEC For the same target accuracy P, the HP75 uses more samples: with P=1E6, the HP71 uses only 127 samples, and the HP75 511 samples. Comparing with the Excel simulations from Albert for the first two results: Code: Intervals Richardson Extrapolations So it turns out that the HP71 and HP75 INTEGRAL are based on the same algorithm but use different criteria for the choice of the number of samples and returned value. The IBOUND value is still not clear for me, both on the HP71 and HP75. With same number of samples and same extrapolation, why is the IBOUND value different and depends on the usersupplied target accuracy P? For instance in the cases above, IBOUND is exactly divided by 10 when P is changed from 1E4 to 1E5  same samples, same extrapolation, same INTEGRAL result. JF 

02052020, 07:52 PM
(This post was last modified: 02052020 08:46 PM by Albert Chan.)
Post: #9




RE: HP71B Integral Questions
(02052020 06:20 PM)JF Garnier Wrote: The IBOUND value is still not clear for me, both on the HP71 and HP75. With same number of samples and same extrapolation, why is the IBOUND value different and depends on the usersupplied target accuracy P? For instance in the cases above, IBOUND is exactly divided by 10 when P is changed from 1E4 to 1E5  same samples, same extrapolation, same INTEGRAL result. I think P is maximum relative error allowed. Let E =  rough expected integral result  → Maximum absolute error ≈ IBOUND = E × P Then, it just walk across the romberg numbers, using the difference as a guide, pick the correct gap. (02032020 01:58 PM)Albert Chan Wrote: Extending rows of 128 intervals, extrapolate to its limits, we have: Say, E = 62563e6 For P=1E5, 63575 < (IBOUND = 625630) < 40399538 → return 62563050740 For P=1E6, 1007 < (IBOUND = 62563) < 63575 → return 62563050648 This implied we need at least 4 romberg numbers, or minimum of 7 sample points. 

02062020, 08:37 AM
(This post was last modified: 02062020 08:57 AM by JF Garnier.)
Post: #10




RE: HP71B Integral Questions
Yes, this is the basic idea, I believe.
I found in the HP71 and HP75 user manuals a significant difference in the way IBOUND is computed in the two machines. This is a simple description, but still helps: HP71: Convergence is determined using J(k) defined as the kth approximation to the integral of F over the same interval of integration. HP75: Convergence is determined using J(k) defined as the kth approximation to the integral of 10^(int(logF)) over the same interval of integration. In both cases, these estimations are compared to the series of integral extrapolations I(k,j)  I(k,j1) < E.J(k) as you described above, using the usersupplied E value. Computing 10^(int(logF)) may look strange and complicate, but it is not with BCD numbers especially in asm: 10^(int(logF)) just means to set the mantissa to 1 and the sign to +, just keeping the exponent. I believe the benefit is to save memory locations, it may come from the original code on the 34C/15C/41C that were quite memorylimited. From this description, the error estimation is using the weighted samples, before any extrapolation. This explains the IBOUND results and the difference between the HP71 and HP75, for instance: HP71: P=?1e5 IA( 127 ) = 62563050744 , IBOUND=625265.300862 . the weighted sample sum was 62526530086, slightly different from the extrapolated value returned by INTEGRAL. HP75: P=?1e3 IA( 127 ) = 62563050559 , IBOUND=21501195.6588 the weighted sample sum (using all mantissa=1) was 21501195659, significantly smaller than the value returned by INTEGRAL. JF 

02062020, 06:53 PM
Post: #11




RE: HP71B Integral Questions
(02052020 08:43 AM)JF Garnier Wrote:(02032020 11:15 PM)Albert Chan Wrote: There was an error in the manual. x = ½(3uu³) substitution ...Actually, this description originates from a HP15C manual (text originally written by W. Kahan) ... and was simplified in the HP71 Math ROM manual. The 71b and 15c manuals' descriptions both quote from Kahan's article in the August 1980 Hewlett Packard Journal description of the 34c [https://www.hpl.hp.com/hpjournal/pdfs/Is...98008.pdf] A couple of notes on the 3/2*u1/2*u³ substitution. 1) The manuals and article say something like Quote:Besides suppressing resonance, the substitution has two additional benefits. First, no sample need be taken from either endpoint of the interval of integration... If I understand correctly (and please correct me if I'm wrong), it's not this substitution that prevents the end points from being evaluated. The reason the endpoints are not evaluated is because the Romberg variation used is based on rectangle midpoint evaluations rather than something like a trapezoid (trapezium) method which does evaluate endpoints. If you were evaluating at the endpoints of 1 and 1, then after this substitution you would still have the endpoints 1 and 1. (This invariance of the endpoints is part of why this substitution is so useful.) 2) Here's an interesting tidbit. I've never seen this in print, but a few years ago I discovered that the ti89/92/nspire also use this substitution even though the intervals are already nonuniform beforehand. The ti83/84 calculators use the standard 15 point GaussKronrod method which uses nonuniform intervals at the following nodes: Code: 0.00000000000000 However, the ti89/92/nspire numeric integration uses the following nodes which didn't match anything I had ever seen. Code: 0.00000000000000 I couldn't figure out where these ti89 values came from until one day on a whim I plugged the GaussKronrod values into 3/2*u1/2*u³ and, voila!, out came the ti89 values. So then the question was "Why make this substitution if the intervals are already nonuniform?" The hp manuals and article go on to say that this substitution allows you to integrate functions whose slope is infinite at an endpoint. While this is not required since the algorithm is not evaluating at the endpoints anyway, experimentation shows that it indeed gives better results when integrating functions like sqrt(1u) or sqrt(1u^2) from 1 to 1 where the slope is infinite at one or more endpoints. By not using the standard GaussKronrod nodes, the algorithm is not as accurate on higher order polynomials (degree > 14), but is more accurate on functions with square roots. Seems like a good trade off. Why this happens is illustrated by the semicircular graph of √(1X^2) before and after the substitution. 

02062020, 11:16 PM
Post: #12




RE: HP71B Integral Questions
(02062020 06:53 PM)Wes Loewer Wrote: A couple of notes on the 3/2*u1/2*u³ substitution. A simple test showed that HP71B INTEGRAL do extrapolations from trapezoids, not midpoints rule. >10 DISP INTEGRAL(1, 1, 1E5, 1/SQRT(1IVAR^2)), IBOUND >RUN 3.14156045534 3.1415397954E5 this failure to converge (65535 sample points !) is due to missing end points evaluation. \(\large \int _{1} ^1 f(x) dx = \int _{1} ^ 1 {3 \over 2}(1u^2) f \left({ u (3u^2) \over 2} \right) du = \int _{1} ^ 1 g(u) du \) \(f(x) = {1\over\sqrt{1x^2}}\quad → \quad g(u) = {3 (1u^2) \over \sqrt{(4u^2)(1u^2)^2}}\) \(\displaystyle{\lim_{u^2 \to 1^} g(u)} = \displaystyle{\lim_{u^2 \to 1^}{3\over\sqrt{4u^2}}} = \sqrt3 ≠ 0\) Quote:So then the question was "Why make this substitution if the intervals are already nonuniform?" It does matter, even if gaussian quadrature do not evaluate end points. utransformation changes the overall shape of the integrand, affecting all sample points. try above f(x) and g(u) (with 1u^2 cancelled) XCas> f(x) := 1/sqrt(1x*x) XCas> g(x) := 3/sqrt(4x*x) XCas> time(gaussquad(f(x), x = 1 .. 1)) → 0.0065 XCas> time(gaussquad(g(x), x = 1 .. 1)) → 0.00054 

02072020, 03:49 AM
(This post was last modified: 02082020 10:35 AM by Wes Loewer.)
Post: #13




RE: HP71B Integral Questions
(02062020 11:16 PM)Albert Chan Wrote:(02062020 06:53 PM)Wes Loewer Wrote: If I understand correctly (and please correct me if I'm wrong), it's not this substitution that prevents the end points from being evaluated. The reason the endpoints are not evaluated is because the Romberg variation used is based on rectangle midpoint evaluations rather than something like a trapezoid (trapezium) method which does evaluate endpoints. But the standard trapezoid method evaluates the endpoints while the rectangle midpoint method does not. As you said, the endpoints are not evaluated, which points to the midpoint method. The 50g definitely uses Romberg extrapolation of midpoints (see edit below) and gives very similar results as above. Setting the number format to FIX 5, you get: result = 3.14156045554, IERR = 3.14153979546E5 I could be wrong, but I would be very surprised if they changed up the algorithm for the 71b. (02062020 11:16 PM)Albert Chan Wrote:Quote:So then the question was "Why make this substitution if the intervals are already nonuniform?" Right, that's just what I was trying to explain and demonstrate in the sample graph. EDIT: Ack! It seems that my statement "The 50g definitely uses Romberg extrapolation of midpoints" was definitely wrong. :( Mea culpa 

02072020, 08:14 AM
(This post was last modified: 02072020 11:09 PM by Albert Chan.)
Post: #14




RE: HP71B Integral Questions
(02072020 03:49 AM)Wes Loewer Wrote: But the standard trapezoid method evaluates the endpoints while the rectangle midpoint method does not. As you said, the endpoints are not evaluated, which points to the midpoint method. The 50g definitely uses Romberg extrapolation of midpoints and gives very similar results as above. Setting the number format to FIX 5, you get: Sorry for the confusion. I thought 65535 sample points failed convergence is enough to show INTEGRAL uses trapezoids, not midpoint rectangles. Here are the actual numbers, and its extrapolations. Midpoint numbers, for g(u)=3/sqrt(4u^2), u = 1 to 1, converged quickly. Note: table required 127 sample points, since previous iteration numbers cannot be reused. Note: we can get midpoint numbers from trapezoids: M_{n} = 2 T_{2n}  T_{n}, even if T_{1} were wrong. Code: 1 3 Trapezoids, with assumption of T_{1} = 0, failed convergence. Last result = 3.141560455, matching HP71B INTEGRAL result. Code: 1 0 The reason for skipping end points is just an optimization, which work for most cases. If f(x) is finite at the end points, g(u) has endpoints = 0, because of (1u^2) factor. If f(x) is infinite at the end points, but grow slower than the shrinking (1u^2), we still have 0 at the endpoints. For other cases, well, at least it returned a negative IBOUND ... 

02072020, 08:23 AM
(This post was last modified: 02072020 08:29 AM by JF Garnier.)
Post: #15




RE: HP71B Integral Questions
(02072020 03:49 AM)Wes Loewer Wrote: The 50g definitely uses Romberg extrapolation of midpoints. It's easy to check that the HP71 doesn't use the mid points but the trapezoid nodes, by tracing the function sampled points: 10 DEF FNF(X) @ PRINT X; @ FNF=X @ END DEF 20 I=INTEGRAL(1,1,.00001,FNF(IVAR)) >RUN 0 .6875 .6875 .3671875 .3671875 .9140625 .9140625 the sampled points corresponds to the u values of 0, +/0.5, +/0.25, +/0.75. I would be interested to know how the 48/49/50 series work. Like the HP71 math ROM, or with further improvements? (I could try but I'm not a RPL guy :) JF 

02072020, 01:19 PM
Post: #16




RE: HP71B Integral Questions
(02072020 08:23 AM)JF Garnier Wrote:(02072020 03:49 AM)Wes Loewer Wrote: The 50g definitely uses Romberg extrapolation of midpoints. Yes, those are the points that you would expect using midpoints. These are the same points used on the 50g. n=1: the rectangle goes from 1 to 1 with a midpoint at 0 n=2: rectangles go from 1 to 0 and 0 to 1 with midpoints at 0.5, 0.5 n=4: rectangles go from 1 to 0.5, 0.5 to 0, 0 to 0.5, and 0.5 to 1 with midpoints at 0.75, 0.25, 0.25, 0.75 

02072020, 05:08 PM
(This post was last modified: 02082020 02:33 AM by Albert Chan.)
Post: #17




RE: HP71B Integral Questions
(02072020 01:19 PM)Wes Loewer Wrote: Yes, those are the points that you would expect using midpoints. These are the same points used on the 50g. Since trapezoids and midpoints are using the same points, I tried to tease out the algorithm, using f(x) = x*log(1+x) 10 DEF FNF(X)=X*LOGP1(X) 20 DEF FNG(U)=1.5*(1U*U)*FNF(U*(3U*U)/2) 30 M1=2*FNG(0) 40 M2=FNG(.5)+FNG(.5) 50 M4=.5*(FNG(.75)+FNG(.25)+FNG(.25)+FNG(.75)) 60 T1=0 70 T2=(T1+M1)/2 80 T4=(T2+M2)/2 90 T8=(T4+M4)/2 100 DISP "MIDPOINTS=";(M120*M2+64*M4)/45 110 DISP "TRAPEZOID=";(T1+84*T21344*T4+4096*T8)/2835 120 DISP "INTEGRAL =";INTEGRAL(1,1,1,FNF(IVAR)) >RUN MIDPOINTS= 1.02693666365 TRAPEZOID= .978017069767 INTEGRAL7= .97801706977 For this case, T_{1} = 0 is a valid assumption. >FNG(.9999), FNG(+.9999) 5.40429438145E3 2.07933751591E4 For a detail comparison between these 2 methods: https://math.stackexchange.com/questions...int/674350 

02072020, 05:54 PM
(This post was last modified: 02072020 05:57 PM by JF Garnier.)
Post: #18




RE: HP71B Integral Questions
(02072020 01:19 PM)Wes Loewer Wrote:(02072020 08:23 AM)JF Garnier Wrote: the sampled points corresponds to the u values of 0, +/0.5, +/0.25, +/0.75. I understand that the 7 values are a single (combined) set of samples, with 8 regions from 1 to 0.75, 0.75 to 0.5 etc. For these 8 regions, midpoints would be 0.875, 0.625, etc. The key point may be that from from one iteration to another, the HP71 doesn't compute a new integral forgetting the previous samples, but combines them with new samples (according to the HP articles/manuals). JF 

02072020, 08:12 PM
Post: #19




RE: HP71B Integral Questions
(02062020 11:16 PM)Albert Chan Wrote: A simple test showed that HP71B INTEGRAL do extrapolations from trapezoids, not midpoints rule. The 50g also chokes on this integral: FIX 3: result = 3.13952114413, IERR = 3.13821310172E3, N = 1023 FIX 4: result = 3.14133371247, IERR = 3.14116984018E4, N = 8191 FIX 5: bails out after N = 10259 It's very well possible that I have simply believed incorrectly for over a dozen years that the Romberg method hp used involved midpoints rather than trapezoids. I went back and reread the Kahan's original article where he writes, "The simplest method is the midpoint rule, whose nodes all lie in the middles of panels all of the same width" along with the statement that "no sample need be drawn from either end of the interval of integration." I guess I jumped to the conclusion that midpoints were used. I sure would love to see the source code on this. 

02072020, 08:16 PM
Post: #20




RE: HP71B Integral Questions
(02072020 05:54 PM)JF Garnier Wrote: The key point may be that from from one iteration to another, the HP71 doesn't compute a new integral forgetting the previous samples, but combines them with new samples (according to the HP articles/manuals). Good point. I picked up on that when I reread the article. 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)