The Museum of HP Calculators

HP Forum Archive 18

[ Return to Index | Top of Index ]

Solving for an Integral: 28S solution
Message #1 Posted by Marcus von Cube, Germany on 28 Feb 2008, 4:42 a.m.

I'm referring here to Chuck's arc length problem which can be found here: Link to archive.

The 28S solver and integrator work differently than in the later machines. Notably, integration can't be used as an algebraic function because it returns two values to the stack, the solution and the error.

Failed algebraic attempt

My first attempt was to create a wrapper to make integration an algebraic function:

<<
  << -> fn var xl xh eps 
    << fn
       var xl xh 3 ->LIST
       eps 
       S
       DROP
    >>
  >>
>>
'INTFN' STO
The local variables fn, var, xl and xh aren't neccessary if you want to execute INTFN from the stack, because the values are fetched from the stack to fill the variables and are than put back in the same order. But the overhead is required by the algebraic syntax.

INTFN can be called like this:

'INTFN(v/(1+4*SQ(T)),T,0,2.3,.0001)' 
EVAL
Due to the algebraic overhead, the function is relatively slow.

Now enter the solver menu:

'INTFN(v/(1+4*SQ(T)),T,0,X,EPS)-6'
STEQ
SOLVR
The resulting variables screen shows 'EPS', 'T' and 'X' as variables. In the parameterlist of INTFN, 'T' is a formal parameter, but the solver logic fails to see this and complains about an unset variable 'T' when I try to solve for 'X'. If I assign a value to 'T', integration fails bacause T is no longer accepted as a formal variable.

So this is a dead end. :-(

RPL solution

I decided that I had to go the RPL way. Integration works without a formal variable if the integrand is a program. Here is my solution:

<< -> x eps
  <<
    << SQ 4 * 1 + v/ >>
    { 0 x }
    eps
    S
    DROP
  >>
>>
'FN' STO

FN can be called as an algebraic:

'FN(2.3,.0001)'
EVAL
This is faster than the original INTFN version.

Now enter the solver menu:

'FN(X,EPS)-6'
STEQ
SOLVR
The resulting variables screen shows 'EPS' and 'X' as a variables. You can now set the accuracy for the integration and solve for 'X'. The result appears within seconds.

:)

I didn't try the ROOT command but it should work likewise.

Marcus

Edited: 28 Feb 2008, 4:43 a.m.

      
Re: Solving for an Integral: 28S solution
Message #2 Posted by Karl Schneider on 28 Feb 2008, 3:25 p.m.,
in response to message #1 by Marcus von Cube, Germany

Hi, Marcus --

Again, good work! You provided a solution for the HP-28C/S, which has different protocol and lesser capabilities for integration than the successor HP-48 series. I'd mentioned this in a post from 18 Feb 2004 within the thread "RPL challenge -- area between two curves" (Archive #14):

Quote:
Arnaud's program can't even be run on a 28C because its integration accepts arguments in a different format (dummy variable and limits enclosed in a list), which I don't know how to automate. It also didn't accept value-loaded variable names for the limits.

You also demonstrated RPL-based integration of a program instead of an algebraic expression -- something I'd previously wondered how to do. That capability is very important, as some things aren't easy to incorporate into algebraic expressions. On the RPN side, the HP-32SII and its successors offer a choice of defining the integrand function by keystroke program or by algebraic expression; all other RPN models require a keystroke program.

I trust that you noticed my solution for the TI-82 (14 Feb 2008), which I based upon your TI-92/Voyage solution. BTW, I find those one-liners clearer than the HP-71 solution, mainly for the absence of the latter's required variable names "IVAR" and "FVAR". However, the TI-82 (at least) cannot perform nested multiple integrations, while the older HP-71 can.

-- KS

Edited: 28 Feb 2008, 3:48 p.m.

            
Re: Solving for an Integral: Optimization for speed.
Message #3 Posted by Marcus von Cube, Germany on 2 Mar 2008, 10:49 a.m.,
in response to message #2 by Karl Schneider

Karl, thank you for the nice words.

I was recently thinking about how the algorithm could be made a little speedier and came to a solution which works in RPL as well as for the 35s in RPN.

The most time is spend in the numerical integration. The first improvement in this special case is to make use of the fact that the derivative of an integral is simply the integrand. Valentin has already jumped in and has modified his generic Newton solver for the 35s.

Numerical solvers iterate in a way that the X value approximates the solution in ever smaller intervalls. This lead me to another simplification:

Integral from 0 to Xn f(t)dt is the same as the sum of Integral from 0 to Xn-1 f(t)dt and Integral from Xn-1 to Xn f(t)dt.

So, instead of computing the numerical integral between 0 and Xn in each iteration, the program can simply save a copy of the last argument and integral value and add the difference in each step. The integration is faster for smaller intervalls as long as the required accuracy is given as an absolute (not relative) value.

It works out for both the 35s and 28S.

My 35S-Implementation:

      LBL F           CHK=188F, LN=29  "Integrand
      RCL T
      x
      4
      *
      1
      +
      v/
      RTN

LBL I CHK=A2B9, LN=19 "Standard integration FN= F I003: 0 x<>y S FN d T RTN

LBL J CHK=CA57, LN=48 "Modified integration x!=0 GTO J008 STO Y "Initialize for X=0 STO Z FN= F RTN J008: FS? 1 "Flag 1 set: use standard intgration GTO I003 RCL Y "Last iteration for X as lower boundary x<>y STO Y "Save new lower boundary S FN d T STO+ Z "Summation RCL Z RTN

LBL S CHK=1954, LN=80 "Solver 0 XEQ J001 "Initialize summation INPUT E "Absolute error in X for solver INPUT X "First guess S006: RCL X STO T "Argument to integrand XEQ F001 "Derivative of integral STO A "Intermediate storage RCL X "Current guess XEQ J001 "Integration, upper boundary on stack 6 - "Subtract right hand side of equation ENTER x<> A "Save result and get derivative 1/x * STO- X "Xn+1 := n - f(n+1) / f'(n+1) ABS RCL E "Maximum absolute error in X x<y? GTO S006 "Next iteration RCL A RCL X "Display value and X STOP GTO S001 "R/S restarts program

Flag 1 decides whether to use the standard (set) or improved (clear) version.

Here is my modfied 28S solution:

<< -> x eps
  << x
     IF
     THEN
       << SQ 4 * 1 + v/ >>
       { XL x } eps
       S
       DROP
       RS +
     ELSE
       0
     END
     x 'XL' STO
     DUP 'RS' STO
  >>
>>
'FN2' STO

If you solve 'FN2(X,EPS)-6' you'll see X, EPS, RS and XL on the menu. RS is the last ReSult, XL ist the Last X value. You need to initialize these by storing 0 to both or storing 0 to X and evaluating the expression once.

The speed improvements are noticable: Depending on the given starting value and accuracy (by display format or variable EPS) the running times can be as little as 25% of the original version.

Edited: 2 Mar 2008, 4:54 p.m.


[ Return to Index | Top of Index ]

Go back to the main exhibit hall