Solving for an Integral: 28S solution



#2

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.


#3

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.


#4

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

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.


Possibly Related Threads...
Thread Author Replies Views Last Post
  [HP-PIRME] BUG Pretty Print, Solution: abs => |...|, ABS => ||...|| CompSystems 2 390 12-13-2013, 09:36 AM
Last Post: CompSystems
  [HP Prime] Tips for Solving Differential Equations More Efficiently Chris Pem10 8 509 11-21-2013, 08:25 PM
Last Post: Chris Pem10
  HP-PRIME CAS SOLVING fabrice48 8 520 10-19-2013, 01:21 PM
Last Post: Han
  HP Prime Solving Nonlinear System of Equations for Complex Results Helge Gabert 11 776 09-30-2013, 03:44 AM
Last Post: From Hong Kong
  Integral Richard Berler 32 1,459 09-22-2013, 08:59 PM
Last Post: Les Koller
  [HP-Prime xCAS] Arrays: Programming Commands (solution to the ambiguity) CompSystems 6 463 08-30-2013, 06:32 PM
Last Post: CompSystems
  HP-Prime: Solution to ambiguity with [ ... ] & another BUGs CompSystems 12 624 08-28-2013, 08:48 AM
Last Post: CompSystems
  Advanced User's Manual and solution to the ambiguity of data types [HP-Prime xCAS] CompSystems 15 973 08-20-2013, 03:37 PM
Last Post: Thomas Klemm
  Riemann's Zeta Function (HP-28S) Gerson W. Barbosa 8 472 02-03-2013, 03:23 PM
Last Post: Gerson W. Barbosa
  [WP34s et al.] Solving the TVM equation for the interest rate Dieter 24 1,202 12-01-2012, 05:53 AM
Last Post: Paul Dale

Forum Jump: