[WP34s] Romberg Integration



#11

Just wanted to advise all fans of the calc that I did a quick and dirty port of the PPC ROM routine IG, and it works beautifully. At this stage I have simply entered it step by step, using where necessary wp34s equivalents (e.g., the 41-series XEQ IND 10 becomes [alpha]XEQ 10). I have yet to streamline it with things like recall register arithmetic and the X~~Y command for comparing rounded interim results. However, it is crazy fast, and with the tough integral discussed by Kahan in his 1980 HP Digest paper, even in ALL 00 the algorithm converges in about a minute or two to a solution that is correct to 14 digits within 1 ULP. In the words of Borat Sagdiyev, high five!

If I am permitted, I will post the entire routine and instructions, but only after I have made enough of my modifications so that it differs enough from the original PPC routine to be openly shared.

Les


#12

Les, write an article in this forum!


#13

I am respectfully boggled that in the months the wp34s has been around no one has mentioned this before. This is not a major task to do a straight port. The wp34s has all of the necessary HP41-style functions and then some. In an initial transfer I had to do next to nothing.

Will look at this more closely as I get to know the wp34s better. I am sure there are efficiencies and whatnot I am missing.

Les


#14

Porting the IG routine has been on my list of things that should be done eventually -- with aim to making this the built in integration routine.

There are definitely efficiencies available and there are other things that would be nice e.g. using the v3 local registers (they can even be resized during the computation if necessary). Then there are the extra fluff required for the built in command -- proper stack on exit and handling infinite values and not numbers.

So while a straight port is very simple, making it integral [sic] to the calculator isn't quite so easy.


- Pauli

#15

Or better, modify it so it is suitable for XROM and we'll include it :-)

ALL nn modes all round to 12 digits. The nn is for display only not for rounding. FIX/SCI and ENG, however, do use their argument for rounding purposes.


- Pauli


#16

Ah, that explains why SCI 11 and ALL 00 each take as long and give the exact same 16 digit final result.


#17

Just out of curiosity, what result do you get when you integrate fractional_portion(x) from 0.0 to 6.4? The same from 0.0 to 6.5? (Sorry, I don't know if the WP34s calls it "FRAC" or "FP" or whatever.)


#18

The integration scheme doesn't like that function at all. Gives 1.28 and 1.625 respectively. Not even close. According to Mathematica, the correct results are 3.08 and 3.125, which is easily seen by inspection (6*1/2 + .4*.4/2 in the first case, and you get the analogy in the second case).

The integral is periodic with discontinuities. The IG Romberg scheme doesn't even seem to get past the first estimate before exiting.

The August 1980 Kahan paper and the advanced integration chapters in the 42s, 15C, 34C, and 35s do at least touch on the "know thy function" theme. In this case, what it easily computed by visual inspection and by hand boggles the calculators. I am sensing that the underlying trapezoidal rule is not fooled by the initial nonlinear distribution of sample points, and the routine gets too successive identical results off the bat and exits. I should note that the algorithm gives exact results quickly within the continuous segments, but not across the discontinuities.

Les


#19

Quote:
The August 1980 Kahan paper and the advanced integration chapters in the 42s, 15C, 34C, and 35s do at least touch on the "know thy function" theme.
If this is what you need: that's easily achieved by quoting the respective texts ;-)
Quote:
In this case, what it easily computed by visual inspection and by hand boggles the calculators.
To me, this problem sounds a bit academic. As known from the vintage calculators you mentioned, discontinuities are a problem. We can't catch every strange case (I doubt any algorithm can), so maybe I have to insert a warning statement in the manual of the kind "Don't use this oven for drying pets!".
#20

At least your results imply that you have a correct port of the IG routine. To the best of my knowledge, this particular example trips up every HP calculator with a built-in integrator. Because of the nonlinear transform, the first three estimates match and the routine refuses to go any further. In a textbook Romberg integration without the transform, the third estimate differs.

The nonlinear transform also bunches points near the endpoints of the interval. This works great when integrating something with a singularity at or near an endpoint, like sqrt(x) from 0 to 1. However, that backfires when the singularity is in the interior of the interval, like sqrt(abs(x)) from -1 to 1.


Possibly Related Threads...
Thread Author Replies Views Last Post
  Integration question and "RPN" mode comment Craig Thomas 16 931 12-05-2013, 02:32 AM
Last Post: Nick_S
  WP34s integration trapped in infinite loop Bernd Grubert 25 1,017 10-17-2013, 08:50 AM
Last Post: Dieter
  HP Prime integration Richard Berler 1 192 10-06-2013, 10:52 PM
Last Post: Helge Gabert
  integration on 39gII emulator Wes Loewer 29 1,200 06-07-2013, 05:58 PM
Last Post: Chris Smith
  WP-34S Integration Richard Berler 15 684 03-08-2013, 02:29 AM
Last Post: Walter B
  HP 34S integration Richard Berler 16 651 02-18-2013, 04:42 PM
Last Post: Marcus von Cube, Germany
  [WP34S] WP34S firmware on the AT91SAM7L-STK dev kit? jerome ibanes 1 203 10-04-2012, 04:59 PM
Last Post: Paul Dale
  [WP34S] Speeding up the Romberg Integration Les Wright 14 571 05-31-2012, 03:39 PM
Last Post: Marcus von Cube, Germany
  New variant for the Romberg Integration Method Namir 8 390 04-18-2012, 07:47 AM
Last Post: Nick_S
  Romberg Integration for 33s, 35s Matt Agajanian 9 373 03-26-2012, 10:00 AM
Last Post: Nick_S

Forum Jump: