Solving for an Integral: 28S solution

[Marcus von Cube, Germany]

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.

[Karl Schneider]

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.

[Marcus von Cube, Germany]

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 X_n f(t)dt is the same as the sum of Integral from 0 to X_n-1 f(t)dt and Integral from X_n-1 to X_n f(t)dt.

So, instead of computing the numerical integral between 0 and X_n 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.