[VA] SRC#001  Spiky Integral

07222018, 12:33 AM
Post: #41




RE: [VA] SRC#001  Spiky Integral
.
Hi, Albert: (07212018 04:00 PM)Albert Chan Wrote: From Valentin Albillo last post, my estimate match correct values, to 15 digits. Comparing values in binary form, accuracy is even better, mine is just 6 ULP over ! You're welcome. Thanks to you for your fresh and effective insight, your contribution to this subject has been truly outstanding and I for one appreciate it. Also, I'm glad if my numeric quadrature results were useful for you to check your remarkably accurate approximations. Quote:You mentioned the problems of it computing I(20000). Sure it does. These are the results my quadrature program outputs for the various ranges you mention: I(20000, 0..Pi/40000) 0 0.000000767448975625012059728677812789736739179600994 1 0.000000767448983276014084555342766291500235218569303 2 0.000000767448983276034720208510846724248989670327563 3 0.000000767448983276034720208518559793576182713599610 4 0.000000767448983276034720208518559793576182713598529 (all 51 decimals digits shown are correct) I(20000, 0..Pi/400000) 0 0.00000076744898327603472020851855979357618264054999631 35064117112150402698116508077139817678384714549561511 1 0.00000076744898327603472020851855979357618264054999631 35038511379081941548835069200608164307257330429915401 2 0.00000076744898327603472020851855979357618264054999631 35038511379081941548835069200608164289806491475047963 (all 106 digits shown are correct) However, I concur with you that the fact that these results are accurate to 51 and 106 decimal digits, respectively, for both shrunk ranges does not necessarily mean that after multiplying them by the correct factor the resulting approximate values for the fullrange (0..2*Pi) integral will be accurate to that many digits, in fact it's quite conceivable that they are not because both results above match only up to 42 digits, namely: 0.000000767448983276034720208518559793576182 Last but not least, using your remarkable shrinking technique, the value on the full integral for N=20,000 indeed can be computed using an HP calc using just its native 12digit precision (!!), as these 3 lines of code I wrote for the HP71B convincingly demonstrates: 1 DESTROY ALL @ SFLAG 1 @ DIM P,I @ DISP 4*INTEGRAL(0,PI/40000,1E8,FNF(IVAR)) 2 DEF FNF(X) @ P=1 @ FOR I=1 TO 20000 @ P=P*COS(I*X) @ IF NOT P THEN END 3 NEXT I @ FNF=P >RUN 3.06979593356E6 (75 min. in Emu71) wich is accurate to 10 digits. Note: That simplest of tricks in the function's definition ("IF NOT P THEN END") halves the computing time. Regards. V. . 

07222018, 04:50 AM
(This post was last modified: 07222018 10:37 AM by Albert Chan.)
Post: #42




RE: [VA] SRC#001  Spiky Integral
(07222018 12:33 AM)Valentin Albillo Wrote: both results above match only up to 42 digits, namely: Hi Valentin, Thanks for the check. The fact that two integral only match 42 digits does not mean anything. The spike formula asked for full spike area. I chopped the spike only for speedier computation, trading accuracy for time. Plotting N vs spike formula digits accuracy, and it followed a straight line, already reaching 15 digits accuracy for N = 60 You might like to try the spike formula again, for N = 1000. For I(1000), you have a correct value to match against ... 

07222018, 02:01 PM
Post: #43




RE: [VA] SRC#001  Spiky Integral
Using Python + numpy, I managed to calculate I(1000)
Code: import numpy >>> spike.spike(1000) 46770858699105378013047692849647150249048020026391352159837485075187255449381044 44575512200800626116981293400844929213580025428715543816767352635000639866305236 49271451391806356678763700553308953903873511563867155215734010997332808966175771 715884247295296277348179194597363883854664431808932677416L Above calculation take 49 seconds on my laptop I(1000) = b[0] / 2^1000 * (2 Pi) = 0.000274258153608378926807734432669808007979394749673091726358234027755841714506 72423455696454538012082538178315765975675889323840397403322977190964502744011004 81739611552026042903356881417709016110968635764441594831973619603002437175558485 42910006760212326308258399120935101281203608176357009606892564575924775067719 ... Rounded to 45 digits, it match Valentin value exactly Is it possible to explain how the integration code work ? How does it handle so many spikes ? 

07222018, 09:06 PM
(This post was last modified: 07222018 09:07 PM by Valentin Albillo.)
Post: #44




RE: [VA] SRC#001  Spiky Integral
,
Hi, Albert: Albert Chan Wrote:Hi Valentin, Thanks for the check. You're welcome. Quote:The fact that two integral only match 42 digits does not mean anything. The spike formula asked for full spike area. [...]You might like to try the spike formula again, for N = 1000. For I(1000), you have a correct value to match against ... I'll do it soon, when I get home in a few hours. Quote:Rounded to 45 digits, it match Valentin value exactly Of course. That's why I stated "(all 48 decimal digits shown are correct)" ... :) ... (the 48 decimal digits mentioned include the three initial zeros as well). Quote:Is it possible to explain how the integration code work ? How does it handle so many spikes ? I developed the code myself as a program for the HP71B, good for 12 digits and easily surpassing the speed of the assemblylanguage Math ROM's INTEGRAL keyword, then ported it nearly verbatim to a multiprecision environment, which is the version I've used to produce the results I've posted in this thread. The code was intended to be included in one of my PDF Datafile articles some 10 years ago, as explaining it and providing relevant examples would take at least some 1012 pages. Alas, for reasons which I won't repeat here the article was never submitted for publication and remains unpublished to this day. If (and when) I can find some adequate place where to publish it, I'll post an announcement here for interested people (like yourself) to download it for free. Thanks for your interest and regards. V. . 

07242018, 12:13 AM
(This post was last modified: 03272019 06:02 PM by Albert Chan.)
Post: #45




RE: [VA] SRC#001  Spiky Integral
Using Gaussian Quadrature for spike formula = 4*I(1000, 0 to Pi/2000), against exact I(1000):


« Next Oldest  Next Newest »

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