41C/CV root finders

[Dave Britten]

I threw together a TVM program - mostly a conversion of the 32S one I did recently, with some local alpha label and flag 22 tricks to make it imitate a 12C - and I'm at the point where I need a numeric solver to compute i%.

I've played around with integrating two different solver routines:

Equation Solver for the HP-41C/CV/CX
Root Finder for the HP-41C/CV/CX (From the Standard Applications module.)

Both of them seem to deteriorate in various ways if the initial guesses aren't very good, given how rapidly the TVM equation starts to grow as i% increases, along with regions where certain subexpressions can be undefined. The first one seemed not too bad if supplied with fairly accurate initial guesses.

Does anybody have any recommendations for good root finders, or techniques for obtaining initial guesses for i%? I know the Advantage module has a general purpose solver that's based on the 15C, but it's also got a TVM solver, which would eliminate the point in writing my own anyway.

[Thomas Klemm]

(05-20-2015 07:07 PM)Dave Britten Wrote:  Does anybody have any recommendations for good root finders, or techniques for obtaining initial guesses for i%?

This post in a recent thread might give you the answer.

Cheers
Thomas

[Dieter]

(05-20-2015 07:07 PM)Dave Britten Wrote:  Does anybody have any recommendations for good root finders, or techniques for obtaining initial guesses for i%?

Welcome back to the wonderful world of TVM solvers. ;-) This subject has been discussed here numerous times, including some basic facts that most probably will not change: there is no simple way of getting a decent guess for all possible combinations of the four other TVM variables. Their signs, and even their values, influence the way a first guess for the interest rate may be determined.

Last time this subject has been discussed in December 2014: have look at this thread, especially post #7 and #14 ff.

That thread also includes a link to the old forum where the problems of solving for the interest rate has been discussed in detail, even with graphs showing different scenarios and their special problems. Here is a link.

I fear things haven't changed much since then. ;-)

Dieter

[Dave Britten]

(05-20-2015 07:36 PM)Thomas Klemm Wrote:  

(05-20-2015 07:07 PM)Dave Britten Wrote:  Does anybody have any recommendations for good root finders, or techniques for obtaining initial guesses for i%?

This post in a recent thread might give you the answer.

Cheers
Thomas

Thanks Thomas, that linked journal article looks useful. I'll pore over it more closely this evening. Whatever they were doing on the HP-80 must have been designed specifically for a slow, space-constrained machine, so that would be a great start.

(05-20-2015 07:44 PM)Dieter Wrote:  

(05-20-2015 07:07 PM)Dave Britten Wrote:  Does anybody have any recommendations for good root finders, or techniques for obtaining initial guesses for i%?

Welcome back to the wonderful world of TVM solvers. ;-) This subject has been discussed here numerous times, including some basic facts that most probably will not change: there is no simple way of getting a decent guess for all possible combinations of the four other TVM variables. Their signs, and even their values, influence the way a first guess for the interest rate may be determined.

Last time this subject has been discussed in December 2014: have look at this thread, especially post #7 and #14 ff.

That thread also includes a link to the old forum where the problems of solving for the interest rate has been discussed in detail, even with graphs showing different scenarios and their special problems. Here is a link.

I fear things haven't changed much since then. ;-)

Dieter

Yeah, I think I was reading some of those threads last night while playing around with the calculation. It's the e^(ln(1+i%)*n) parts that seem to grow rapidly, so I tried solving just those for i%, coming up with e^(x/n)-1=i%, then used x=1 and x=10 to get starting values for i%. The "Susan problem" converged very rapidly, but solving for i% in my own mortgage caused it to oscillate wildly around the real interest rate over a range of about two orders of magnitude. It converged eventually, but took an unusually long time.

But I think both of the solvers I tried were based on the secant method. Maybe Newton-Raphson performs better (its usage in the HP-80 seems to suggest so).

On a side note, it's rather strange that the 41C has exotic functionality like ln1+X, e^x-1, mod, and extensive alpha capability, but lacks a root finder, present on earlier models like the 34C. Not to mention the lack of combinatorics and hyperbolic trig (good thing it already has e^x-1, since I don't know a good way to calculate it without sinh).

[Massimo Gnerucci]

(05-20-2015 08:08 PM)Dave Britten Wrote:  On a side note, it's rather strange that the 41C has exotic functionality like ln1+X, e^x-1, mod, and extensive alpha capability, but lacks a root finder, present on earlier models like the 34C. Not to mention the lack of combinatorics and hyperbolic trig (good thing it already has e^x-1, since I don't know a good way to calculate it without sinh).

The nice thing with the 41 is that you don't need to have everything on board: if you want something exotic just key in the program from a Users' Library Solutions book or, better, plug in the required module. ;)

[Dave Britten]

(05-20-2015 08:40 PM)Massimo Gnerucci Wrote:  

(05-20-2015 08:08 PM)Dave Britten Wrote:  On a side note, it's rather strange that the 41C has exotic functionality like ln1+X, e^x-1, mod, and extensive alpha capability, but lacks a root finder, present on earlier models like the 34C. Not to mention the lack of combinatorics and hyperbolic trig (good thing it already has e^x-1, since I don't know a good way to calculate it without sinh).

The nice thing with the 41 is that you don't need to have everything on board: if you want something exotic just key in the program from a Users' Library Solutions book or, better, plug in the required module.

Yup, that's the big advantage. If it didn't have the extensive memory and expandability, I probably wouldn't give it a second thought, and just stay with my 32S. It's still odd that some of those basics were omitted, though. Linear regression is missing as well. I wonder if they had to keep ROM size down to allow for enough RAM at a reasonable price point.

[Didier Lachieze]

(05-20-2015 08:08 PM)Dave Britten Wrote:  On a side note, it's rather strange that the 41C has exotic functionality like ln1+X, e^x-1, mod, and extensive alpha capability, but lacks a root finder, present on earlier models like the 34C.

Well, while being the last HP calculator with red LED display, the 34C was introduced after the 41C, not before: it was introduced along with the 33C and the 38C on August 15, 1979, one month after the 41C introduction on July 16, 1979 (source: PPC Journal V6N4 & V6N5).

[Dieter]

(05-20-2015 07:07 PM)Dave Britten Wrote:  I've played around with integrating two different solver routines:

Equation Solver for the HP-41C/CV/CX
Root Finder for the HP-41C/CV/CX (From the Standard Applications module.)

Both of them seem to deteriorate in various ways if the initial guesses aren't very good...

FTR: both programs use a regula falsi method. The first one uses the simple standard method, the second one seems to implement a modified version, the Illinois method. The latter should converge to a solution within an interval if (!) the two initial guesses bracket the solution, i.e. f(x1) and f(x2) have opposite signs. Further modifications can be found in the German Wikipedia article.

(05-20-2015 07:07 PM)Dave Britten Wrote:  ...given how rapidly the TVM equation starts to grow as i% increases, along with regions where certain subexpressions can be undefined.

That's what I already said in earlier discussions: there are various ways to write the TVM equation. Some are more suitable for solving, others less so. This point has even been discussed in the HP Journal article on the then-new HP92 (cf. Oct 1977 issue, p. 24/25). Here the Newton method was chosen, along with an initial guess that is claimed to be "correct to five places over all ranges of PV, FV, PMT, and i, and with n as large as 10⁸". Which sounds truly amazing.

Dieter

[Dave Britten]

I knocked out a simple lazy Newton's method solver (lazy in that it approximates the 1st derivative using f(1.000001x) and f(.999999x)), and it seems functional but definitely suboptimal. Here's the code for it, if you want to have a chuckle at my background obviously not being in numeric methods.

Code:


LBL "FN="
"FN="
AON
PROMPT
ASTO 20
AOFF
RTN

LBL "SOLVE"
"REG="
AOFF
PROMPT
STO 21
LBL "PSOLVE"
RCL IND 21
STO 22
LBL 00
CLA
ARCL X
AVIEW
STO IND 21
XEQ IND 20
STO 23
X=0?
GTO 04
RCL 22
1
1 E-6
-
*
STO IND 21
XEQ IND 20
STO 24
RCL 22
1
1 E-6
+
*
STO IND 21
XEQ IND 20
STO 25
RCL 24
-
X=0?
GTO 09
RCL 22
2 E-6
*
X<>Y
/
RCL 23
*
CHS
RCL 22
+
STO 22
GTO 00
LBL 04
RCL 22
RTN
LBL 09
RCL 23
RCL 22
RTN
END

LBL 09 is where it would check for extrema or discontinuities if I had implemented that part yet. I fed it numbers similar to my mortgage loan, and it obtained i%, but it wasn't exactly fast. Evaluating the function three times for every iteration really takes its toll on speed.

I'll have to study the HP-80 article more so as to have a better approximation of its algorithm. There's also a 19C/29C i% solver that uses Newton's method that's on the museum DVD (part of the general application book, I think), and I might be able to use some of that. I don't even necessarily need a general purpose solver, just one that's suited for handling the TVM equation.

[Dieter]

(05-21-2015 12:51 PM)Dave Britten Wrote:  I knocked out a simple lazy Newton's method solver (lazy in that it approximates the 1st derivative using f(1.000001x) and f(.999999x)), and it seems functional but definitely suboptimal.

Well.... ;-)

(05-21-2015 12:51 PM)Dave Britten Wrote:  I'll have to study the HP-80 article more so as to have a better approximation of its algorithm.

I would not recommend the HP80 as a reference. It used a different TVM paradigm with unsigned variables that was changed with the HP92, which essentially used the method that was later implemented in many classic HP financial calculators.

(05-21-2015 12:51 PM)Dave Britten Wrote:  I don't even necessarily need a general purpose solver, just one that's suited for handling the TVM equation.

I already suggested the method in post #7 of the older thread linked in my previous post. Here is a quick and dirty implementation of that idea. Try it and see what you get.

Code:

Registers:
R01 = n
R02 = i%
R03 = PV
R04 = PMT
R05 = FV

R07 = i%/100
R08 = used for variable a = (1+i%/100)^(-n)
R09 = used for variable b = pmt / (i%/100) * (1 - a)

Flag 01 clr = END mode
Flag 01 set = BEGIN mode

HP-41 program:

01 LBL"I%"
02 RCL 01
03 RCL 04
04 *
05 RCL 03
06 +
07 RCL 05
08 +
09 X=0?      // in this case we're already done
10 GTO 99
11 1
12 FC? 01
13 CHS
14 RCL 01
15 +
16 RCL 05
17 *
18 1
19 FC? 01
20 CHS
21 RCL 01
22 -
23 RCL 03
24 *
25 +
26 /
27 X>0?      // can be omitted, but seems to improve the guess
28 ST+ X
29 STO 07    // save first guess as i%/100
30 LBL 01    // start iteration
31 RCL 07
32 LN1+X
33 RCL 01
34 *
35 CHS
36 E^X
37 STO 08    // save (1+i%/100)^(-n)
38 LASTX
39 E^X-1     // evaluate (1+i%/100)^(-n) - 1 accurately
40 CHS
41 RCL 04
42 *
43 RCL 07
44 /
45 STO 09
46 RCL 09
47 RCL 07
48 *
49 FC? 01
50 CLX
51 +
52 RCL 03
53 +
54 RCL 08
55 RCL 05
56 *
57 +      // = TVM(i)
58 RCL 07
59 1/x
60 1
61 FC? 01
62 CLX
63 +
64 RCL 04
65 *
66 RCL 05
67 -
68 RCL 01
69 *
70 RCL 08
71 *
72 RCL 07
73 1
74 +
75 /
76 RCL 09
77 RCL 07
78 /
79 -         // = TVM'(i)
80 /         // = correction term for i%/100
81 ST- 07    // adjust i%/100
82 RCL 07
83 /
84 ABS       // = abs(relative error)
85 100
86 +
87 LASTX
88 X#Y?      // is relative error not yet below 5 E-8 ?
89 GTO 01    // then do another iteration
90 RCL 07
91 *         // calculate i% (here Y still holds the 100 from step 85/87)
92 LBL 99
93 STO 02    // save result in R02
94 END       // and exit

This is essentially untested code, so beware. ;-) But it seems to work.
Try the standard example:

  10 STO 01 (n)
  50 STO 03 (PV)
–30 STO 04 (PMT)
400 STO 05 (FV)

XEQ"I%" => 14,4359

The program evaluates a first guess of 9,8361% and finally returns 14,4359% after five iterations. With flag 01 set (BEGIN) the result is 10,2051%.

Dieter

[Dieter]

(05-21-2015 07:36 PM)Dieter Wrote:  I already suggested the method in post #7 of the older thread linked in my previous post. Here is a quick and dirty implementation of that idea. Try it and see what you get.

FTR, here is another version with a different definition of the TVM equation. It is shorter, runs faster and may require slightly less iterations here and there. Maybe you want to try this as well. Lines 01...29 are the same as before, here is the new code from line 30:

Code:

.. ...
30 LBL 01
31 RCL 04
32 RCL 03
33 FS? 01
34 +
35 RCL 04
36 RCL 07
37 /
38 +
39 STO 08    // a = pv + pmt/(i%/100) ( + pmt for END)
40 RCL 07
41 LN1+X
42 RCL 01
43 *
44 E^X-1
45 STO 09    // b = (1+i%/100)^n - 1
46 *
47 RCL 03
48 +
49 RCL 05
50 +         // = TVM
51 RCL 09
52 1
53 +
54 RCL 07
55 1
56 +
57 /
58 RCL 01
59 *
60 RCL 08
61 *
62 RCL 04
63 RCL 09
64 *
65 RCL 07
66 x^2
67 /
68 -         // = TVM'
69 /         // = correction term
70 ST- 07    // adjust i%/100
71 RCL 07
72 /
73 ABS       // = abs(relative error)
74 100
75 +
76 LASTX
77 X#Y?      // is relative error not yet below 5 E-8 ?
78 GTO 01    // then do another iteration
79 RCL 07
80 *         // calculate i% (here Y still holds the 100 from step 85/87)
81 LBL 99
82 STO 02    // save result in R02
83 END       // and exit

The standard problem was solved with just four iterations (instead of five before) in about 10 seconds.

Dieter

[Ángel Martin]

(05-20-2015 07:07 PM)Dave Britten Wrote:  Does anybody have any recommendations for good root finders, or techniques for obtaining initial guesses for i%?

Here's the formula I used in the TVM Module for the first derivative (as it uses Newton's method):

f ’ (i) = (PMT / i^2 ) * [ (1+i)^(-n) - 1 ] + n * [PMT (1 + ip)/i – FV ] * (1+i)^-(n+1)

with initial value:
i0 = [ abs( PV + n*PMT + FV) ]^1/n

Seems to hold up quite well but again this is a MCODE implementation using the internal 13-digit precision routines from the OS.

Note: edited to add the missing exponent signs.


Cheers,
AM

[Dieter]

(05-22-2015 01:14 PM)Ángel Martin Wrote:  Here's the formula I used in the TVM Module for the first derivative (as it uses Newton's method):

f ’ (i) = (PMT / i2 ) * [ (1+i)-n - 1 ] + n * [PMT (1 + ip)/i – FV ] * (1+i)-(n+1)

with initial value:
i0 = [ abs( PV + n*PMT + FV) ]1/n

Hmmm... looks like there are some characters missing. ;-)
How did you arrive at the formula for the inital guess?

For the record, here are the formulas I used for the posted 41C programs.
In all cases, k=0 (END) resp. 1 (BEGIN).

The initial guess is calculated as follows:

Code:

s = 2*k - 1
i = 2 * (pv + n*pmt + fv) / (fv * (s+n) + pv * (s-n))

That's essentially the value of the next Newton iteration if you start at i=0. ;-)

The first program uses a TVM formula close to the definition in the 12C manual. It is evaluated this way:

Code:

q = 1 + i
a = q^(-n)
b = pmt / i * (1 - a)

tvm = pv + b * (1 + k*i) + fv * a
tvm_deriv = n * a / q * (pmt * (1/i + k) - fv) - b/i

The second, improved version uses the following formulas:

Code:

q = 1 + i
a = pv + pmt * k + pmt / i
b = q^n - 1

tvm = a * b + pv + fv
tvm_deriv = a * n * q^(n-1) - pmt * b / i^2
          = a * n * (b+1)/q - pmt * b / i^2

This way slow operations like powers and logs are reduced to a minimum.
Of course the actual implementations use e^x–1 and ln(1+x) for better accuracy.

Addendum:
Here's another way of calculating a first guess for the interest rate. Again, it's simply the first approximation using Newton's method starting at i=0. But this time based on the second TVM definition:

Code:

s = 2*k - 1
i = 2 / n * (pv + n*pmt + fv) / (pmt * (s+n) + 2 * pv)

For the test case this leads to an initial estimate of 17,647% (END) resp. 13,044% (BEGIN) while the final solutions are 14,436% resp. 10,205%.

Dieter

[Dieter]

(05-22-2015 01:14 PM)Ángel Martin Wrote:  Seems to hold up quite well but again thit is a MCODE implementation using the internal 13-digit precision routines from the OS.

Ángel, how does your MCODE version behave if there is no solution? Does it loop forever, did you provide an iteration counter, or what will happen? And – please excuse my dumb question – is MCODE interruptable by the user?

You may try the standard example (n=10, PV=50, PMT=–30, FV=400) with n=8 or less. In this case no solution exists.

Some time ago I tried to handle such cases in an extensive 34s program that was able to detect whether a solution exists or not (well, at least in most cases ;-)). This included finding a local extreme of the TVM function and checking for sign changes, so that even multiple solutions could be found. But covering at least the most common cases is not trivial at all, cf. earlier discussions in the old forum.

Dieter

[Ángel Martin]

(05-24-2015 07:10 PM)Dieter Wrote:  Ángel, how does your MCODE version behave if there is no solution? Does it loop forever, did you provide an iteration counter, or what will happen? And – please excuse my dumb question – is MCODE interruptable by the user?

You may try the standard example (n=10, PV=50, PMT=–30, FV=400) with n=8 or less. In this case no solution exists.

Yes it is interruptible (the code checks for a key press event) but no, there is no counter due to available space restriction (it is also in the SandMath, which was very crowded already). So the example above just keeps iterating until the user stops it. Not perfect, but then... who said "the better is the enemy of the good"?

As to whether there are missing terms in the formula for the derivative, well I hope not; barring a sinister transcription error that's the expression used in the code. I'll double check my notes.

Cheers,
'AM

[Dieter]

(05-25-2015 05:35 AM)Ángel Martin Wrote:  Yes it is interruptible (the code checks for a key press event) but no, there is no counter due to available space restriction (it is also in the SandMath, which was very crowded already). So the example above just keeps iterating until the user stops it. Not perfect, but then... who said "the better is the enemy of the good"?

Voltaire?

(05-25-2015 05:35 AM)Ángel Martin Wrote:  As to whether there are missing terms in the formula for the derivative, well I hope not; barring a sinister transcription error that's the expression used in the code. I'll double check my notes.

I think there are some power symbols "^" missing. Both in the derivative and in the initial guess.

I assume the latter should read [abs(PV + n*PMT + FV)]^1/n. For the sample case in this thread this yields an initial i of 1,65 or 165%. ?!? Also the value can never be negative. That's why I asked how you developed this formula.

BTW if space is an issue there is a very simply solution: simply set the inital i to a value very close to zero, e.g. 10^–4. A value of exactly zero will not work because the standard TVM formulas divide by i, but a value very close to this would return (approximately) the estimate I posted - right after the first iteration. So all you lose is a bit of execution speed – i.e. the difference between the time required for one iteration loop and the calculation of the initial estimate. ;-)

Dieter

[Ángel Martin]

(05-25-2015 07:19 PM)Dieter Wrote:  I think there are some power symbols "^" missing. Both in the derivative and in the initial guess.

You're right - sorry I didn't notice those before, The exponent signs were not handled well in the copy-paste action; I have corrected the expressions in the initial post now.

The initial guess just look intuitively appropriate to me, I confess I didn't derive it from any algorithm or deductive process... but it just looked right - yes I know, math heresy but sometimes one must follow the hunch even in math.

(05-25-2015 07:19 PM)Dieter Wrote:  BTW if space is an issue there is a very simply solution: simply set the inital i to a value very close to zero, e.g. 10^–4. A value of exactly zero will not work because the standard TVM formulas divide by i, but a value very close to this would return (approximately) the estimate I posted - right after the first iteration. So all you lose is a bit of execution speed – i.e. the difference between the time required for one iteration loop and the calculation of the initial estimate. ;-)

I fail to see how reducing the initial guess to a near-zero value will solve the problem when there is not a solution to be found...which is the only problem left. Again, the initial guess I use works rather well for all other *solvable* cases, it's sweet and short to obtain and doesn't require too many iterations to converge.

The SandMath manual has some further details in case it helps.

Thanks Dieter for the suggestions.

'AM

[Dieter]

(05-26-2015 05:06 AM)Ángel Martin Wrote:  The exponent signs were not handled well in the copy-paste action; I have corrected the expressions in the initial post now.

OK. But back to the initial guess: does it really evaluate to 165% for the sample case, or did I understand something wrong?

(05-26-2015 05:06 AM)Ángel Martin Wrote:  I fail to see how reducing the initial guess to a near-zero value will solve the problem when there is not a solution to be found...which is the only problem left.

No, of course an estimate close to zero will not solve cases where no (real) solution exists. But it saves memory space at the cost of merely one more iteration.

(05-26-2015 05:06 AM)Ángel Martin Wrote:  Again, the initial guess I use works rather well for all other *solvable* cases, it's sweet and short to obtain and doesn't require too many iterations to converge.

Please forgive me if I get back to this point once again, but I simply want to know: How many iterations are required for the example posted in this thread until the result (14,436%) is returned? Does the algorithm actually return the lower one of the two solutions (the other one being 53,172%), although the first estimate is as high as 165%? Since |PV+n*PMT+FV| usually is > 1, the estimate will (almost) always be > 100%.

And finally... may I ask this...? Would you mind trying your algorithm with a first estimate of 1E–4? How does it behave?

Dieter

[Ángel Martin]

(05-26-2015 05:39 PM)Dieter Wrote:  OK. But back to the initial guess: does it really evaluate to 165% for the sample case, or did I understand something wrong?

Not sure I follow, I didn't try the example before - but I just ran it an here are the results:

10 = (n)
50 = (PV)
–30 =(PMT)
400 = (FV)

I=14.43587035% in end mode, and
I=10.20509732% in begin mode

you can also run it yourself on V41 of course...

(05-26-2015 05:39 PM)Dieter Wrote:  How many iterations are required for the example posted in this thread until the result (14,436%) is returned? Does the algorithm actually return the lower one of the two solutions (the other one being 53,172%), although the first estimate is as high as 165%? Since |PV+n*PMT+FV| usually is > 1, the estimate will (almost) always be > 100%.

And finally... may I ask this...? Would you mind trying your algorithm with a first estimate of 1E–4? How does it behave?

I'd need to put breakpoints in the code to read the initial guess and the number of iterations, let me work on that and I'll get back to you tomorrow.

Cheers,
'AM

[Dieter]

(05-26-2015 09:26 PM)Ángel Martin Wrote:  Not sure I follow, I didn't try the example before - but I just ran it an here are the results:

10 = (n)
50 = (PV)
–30 =(PMT)
400 = (FV)

I=14.43587035% in end mode, and
I=10.20509732% in begin mode

That's close to the exact results (14,43587133% and 10,20509807%). I still wonder why an estimate > 100% leads to these results instead of the second, larger solutions.

(05-26-2015 09:26 PM)Ángel Martin Wrote:  you can also run it yourself on V41 of course...

I could do so if I was more familiar with V41. ;-)

(05-26-2015 09:26 PM)Ángel Martin Wrote:  I'd need to put breakpoints in the code to read the initial guess and the number of iterations, let me work on that and I'll get back to you tomorrow.

Thank you very much. :-)

Dieter