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 X_{n} f(t)dt is the same as the sum of Integral from 0 to X_{n1} f(t)dt and Integral from X_{n1} 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 35SImplementation:
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 "X_{n+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.
