+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html)
+--- Forum: General Forum (/forum-4.html)
+--- Thread: 41C/CV root finders (/thread-3925.html)
RE: 41C/CV root finders - Dieter - 05-29-2015 06:15 PM
(05-29-2015 08:04 AM)Ángel Martin Wrote: The formulas are posted below, but this is certainly strange - I need to revise the implementation, it might be something goofed up in the way I program them...
That's probably true. See below.
(05-29-2015 08:04 AM)Ángel Martin Wrote: f(i) = PV + (1+ip) PMT [ (1- (1+i)^-n) / i] + FV (1+i)^-n
f ’ (i) = (PMT / i^2 ) * [ (1+i)^(-n) - 1 ] + n * [PMT (1 + ip)/i – FV ] * (1+i)^-(n+1)
where p=1 in Begin mode or 0 in End mode.
The formulas look a bit complicated, and I hope they are implemented efficiently, but anyway: they seem to be correct. So I tried Newton's method based on your formulas, and now look what happens even with a less-than-optimal initial guess like the one you posted. The calculations were done in Excel, so the results may slighty differ, but the essential point should be clear:
Code:
# i f(i) f'(i) f(i)/f'(i) new i
---------------------------------------------------------------------------------------------------------------------------
0 1,6504 75667 020 E-02 1,1512 56813 905 E+02 -1,8922 05978 309 E+03 -6,0842 04505 758 E-02 7,7346 80172 779 E-02
1 7,7346 80172 779 E-02 3,6157 27302 052 E+01 -8,3765 28811 450 E+02 -4,3164 98377 119 E-02 1,2051 17854 990 E-01
2 1,2051 17854 990 E-01 9,0500 96755 600 E+00 -4,5257 59747 816 E+02 -1,9996 85635 095 E-02 1,4050 86418 499 E-01
3 1,4050 86418 499 E-01 1,2439 53188 362 E+00 -3,3307 99490 524 E+02 -3,7346 98506 773 E-03 1,4424 33403 567 E-01
4 1,4424 33403 567 E-01 3,6176 85118 950 E-02 -3,1385 44999 486 E+02 -1,1526 63135 160 E-04 1,4435 86066 702 E-01
5 1,4435 86066 702 E-01 3,3398 41367 733 E-05 -3,1327 51375 488 E+02 -1,0661 04828 447 E-07 1,4435 87132 807 E-01
6 1,4435 87132 807 E-01 2,8876 45678 129 E-11 -3,1327 46020 753 E+02 -9,2176 18214 181 E-14 1,4435 87132 808 E-01With the estimate I suggested things look like this:
Code:
# i f(i) f'(i) f(i)/f'(i) new i
---------------------------------------------------------------------------------------------------------------------------
0 8,5714 28571 429 E-02 2,9538 18144 886 E+01 -7,4603 90801 804 E+02 -3,9593 34334 298 E-02 1,2530 76290 573 E-01
1 1,2530 76290 573 E-01 6,9556 69566 821 E+00 -4,2116 88482 615 E+02 -1,6515 15679 646 E-02 1,4182 27858 537 E-01
2 1,4182 27858 537 E-01 8,1076 20873 942 E-01 -3,2621 35750 476 E+02 -2,4853 72005 981 E-03 1,4430 81578 597 E-01
3 1,4430 81578 597 E-01 1,5844 14941 007 E-02 -3,1352 86065 219 E+02 -5,0534 94029 090 E-05 1,4435 86928 000 E-01
4 1,4435 86928 000 E-01 6,4161 18566 221 E-06 -3,1327 47049 439 E+02 -2,0480 80634 972 E-08 1,4435 87132 808 E-01
5 1,4435 87132 808 E-01 1,0231 81539 495 E-12 -3,1327 46020 748 E+02 -3,2660 85194 006 E-15 1,4435 87132 808 E-01
6 1,4435 87132 808 E-01 0,0000 00000 000 E+00 -3,1327 46020 748 E+02 0,0000 00000 000 E+00 1,4435 87132 808 E-01And here are the results for BEGIN mode with the better estimate:
Code:
# i f(i) f'(i) f(i)/f'(i) new i
---------------------------------------------------------------------------------------------------------------------------
0 8,5714 28571 429 E-02 1,2719 70870 681 E+01 -8,6744 78838 720 E+02 -1,4663 36934 276 E-02 1,0037 76550 570 E-01
1 1,0037 76550 570 E-01 1,1771 83138 188 E+00 -7,1156 23484 535 E+02 -1,6543 64007 239 E-03 1,0203 20190 643 E-01
2 1,0203 20190 643 E-01 1,3189 32492 805 E-02 -6,9567 11221 997 E+02 -1,8959 13817 199 E-05 1,0205 09782 025 E-01
3 1,0205 09782 025 E-01 1,7088 07900 513 E-06 -6,9549 08670 105 E+02 -2,4569 81078 499 E-09 1,0205 09806 594 E-01
4 1,0205 09806 594 E-01 -3,9790 39320 257 E-13 -6,9549 08436 533 E+02 5,7211 95838 259 E-16 1,0205 09806 594 E-01
5 1,0205 09806 594 E-01 0,0000 00000 000 E+00 -6,9549 08436 533 E+02 0,0000 00000 000 E+00 1,0205 09806 594 E-01
6 1,0205 09806 594 E-01 0,0000 00000 000 E+00 -6,9549 08436 533 E+02 0,0000 00000 000 E+00 1,0205 09806 594 E-01So in all cases 4 – 6 iterations will do, depending on the initial estimate.
That's why I wonder why 7 to 13 iterations are required before the solver in the TVM ROM converges. This sounds like an error in the implementation. Maybe you can compare the results of every iteration step with the ones in the tables above.
By the way, here are the results for an initial estimate of 0,0001:
Code:
END mode:
# i f(i) f'(i) f(i)/f'(i) new i
---------------------------------------------------------------------------------------------------------------------------
0 1,0000 00000 000 E-04 1,4976 51539 336 E+02 -2,3469 21996 559 E+03 -6,3813 43485 345 E-02 6,3913 43485 345 E-02
1 6,3913 43485 345 E-02 4,8510 92508 382 E+01 -1,0063 21939 891 E+03 -4,8206 16858 364 E-02 1,1211 96034 371 E-01
2 1,1211 96034 371 E-01 1,3094 15813 711 E+01 -5,1222 53398 655 E+02 -2,5563 27678 077 E-02 1,3768 28802 179 E-01
3 1,3768 28802 179 E-01 2,2064 35164 358 E+00 -3,4822 77391 105 E+02 -6,3361 84388 970 E-03 1,4401 90646 068 E-01
4 1,4401 90646 068 E-01 1,0669 34219 175 E-01 -3,1498 41444 676 E+02 -3,3872 63257 263 E-04 1,4435 77909 326 E-01
5 1,4435 77909 326 E-01 2,8895 04148 129 E-04 -3,1327 92347 864 E+02 -9,2234 14217 349 E-07 1,4435 87132 740 E-01
6 1,4435 87132 740 E-01 2,1364 88319 593 E-09 -3,1327 46021 090 E+02 -6,8198 58058 105 E-12 1,4435 87132 808 E-01
7 1,4435 87132 808 E-01 -1,1368 68377 216 E-13 -3,1327 46020 748 E+02 3,6289 83548 895 E-16 1,4435 87132 808 E-01
BEGIN mode:
# i f(i) f'(i) f(i)/f'(i) new i
---------------------------------------------------------------------------------------------------------------------------
0 1,0000 00000 000 E-04 1,4973 51704 270 E+02 -2,6465 92194 474 E+03 -5,6576 59337 907 E-02 5,6676 59337 907 E-02
1 5,6676 59337 907 E-02 4,3448 06469 682 E+01 -1,2732 67705 318 E+03 -3,4123 27550 236 E-02 9,0799 86888 143 E-02
2 9,0799 86888 143 E-02 8,4553 78795 595 E+00 -8,1017 60264 595 E+02 -1,0436 47123 520 E-02 1,0123 63401 166 E-01
3 1,0123 63401 166 E-01 5,6973 99825 850 E-01 -7,0327 39452 499 E+02 -8,1012 52526 604 E-04 1,0204 64653 693 E-01
4 1,0204 64653 693 E-01 3,1404 39865 149 E-03 -6,9553 37692 589 E+02 -4,5151 50815 028 E-06 1,0205 09805 201 E-01
5 1,0205 09805 201 E-01 9,6905 58044 895 E-08 -6,9549 08449 779 E+02 -1,3933 40849 110 E-10 1,0205 09806 594 E-01
6 1,0205 09806 594 E-01 0,0000 00000 000 E+00 -6,9549 08436 533 E+02 0,0000 00000 000 E+00 1,0205 09806 594 E-01So in fact just one more iteration is required.
(05-29-2015 08:04 AM)Ángel Martin Wrote: PS.- You don't need the Library#4 for the TVM$ module - althought it doesn't hurt to have it always configured on V41 of course. You'll need it for the SandMath, SandMatrix and 41Z if you decide to venture there.
OK, good to know. I now am looking a bit into the functions of the Sandmath ROM.
Dieter
RE: 41C/CV root finders - Ángel Martin - 05-29-2015 06:44 PM
(05-29-2015 06:15 PM)Dieter Wrote:(05-29-2015 08:04 AM)Ángel Martin Wrote: f(i) = PV + (1+ip) PMT [ (1- (1+i)^-n) / i] + FV (1+i)^-n
f ’ (i) = (PMT / i^2 ) * [ (1+i)^(-n) - 1 ] + n * [PMT (1 + ip)/i – FV ] * (1+i)^-(n+1)
where p=1 in Begin mode or 0 in End mode.
The formulas look a bit complicated, and I hope they are implemented efficiently, but anyway: they seem to be correct.
I think they are correct and correctly implemented now after I fixed the bug - see previous post. - but you're welcome to improve the method any time ;-)
(05-29-2015 06:15 PM)Dieter Wrote: So in all cases 4 – 6 iterations will do, depending on the initial estimate.
That's why I wonder why 7 to 13 iterations are required before the solver in the TVM ROM converges. This sounds like an error in the implementation. Maybe you can compare the results of every iteration step with the ones in the tables above.
You may have missed my previous post, or perhaps our posts crossed paths...
it is down to 4 or 5 iterations now after I fixed the bug I found in the algorithm.
Cheers,
'AM
RE: 41C/CV root finders - Dieter - 05-29-2015 07:19 PM
(05-29-2015 06:44 PM)Ángel Martin Wrote: I think they are correct and correctly implemented now after I fixed the bug - see previous post. - but you're welcome to improve the method any time ;-)
I am using the TVM and TVM' formulas posted earlier in this thread (cf. post #11 and #13). The implementation uses two intermediate results a and b that appear both in the function and the derivative, so the slow functions like powers and logs are minimized. This way the iteration runs faster and uses less memory.
Let's see if there is a similar way to optimize the implementation of the formulas you use. ;-)
But, more important:
(05-29-2015 06:44 PM)Ángel Martin Wrote: You may have missed my previous post, or perhaps our posts crossed paths...
Yes, it took me about 20 minutes to compile all the data, so I started editing my post before I could read your update. #-)
(05-29-2015 06:44 PM)Ángel Martin Wrote: ...it is down to 4 or 5 iterations now after I fixed the bug I found in the algorithm.
Great – mission accomplished. Now the community waits for an updated TVM 1F. :-)
Dieter
RE: 41C/CV root finders - Dieter - 05-29-2015 08:00 PM
(05-29-2015 06:44 PM)Ángel Martin Wrote: ...but you're welcome to improve the method any time ;-)
There actually is another potential improvement: PMT=0 currently is not handled separately. At least the "Running..." display appears for several seconds on V41, and "impossible cases" like both PV and FV being positive seem to loop forever.
In a FOCAL program I would implement the direct calculation for PMT=0 this way:
Code:
RCL 03
RCL 05
+
CHS
RCL 03
/
ln1+x
RCL 01
/
e^x-1
100
*Sometimes one forgets the most simple solution. ;-)
Dieter
RE: 41C/CV root finders - Dieter - 05-30-2015 12:37 PM
(05-29-2015 07:19 PM)I Wrote: Let's see if there is a similar way to optimize the implementation of the formulas you use. ;-)
There is a similar method like the one suggested in post #13. The calculation of your f(i) and f'(i) could get implemented this way:
Code:
Let
a = expm1(-n * ln1p(i))
b = PMT * (1/i + p) - FV
Then
-f(i) = a * b - PV - FV
f'(i) = a * PMT/i^2 + b * n * (1+a)/(1+i)Afterwards –f(i) is divided by f'(i) and this quotient is the delta_i that is then added (!) to the previous estimate. Now every iteration requires just one ln1+x and one e^x-1, the formulas are shorter and maybe even more accurate.
Regarding cases with no solutions and infinite looping: I think there is a better way than a simple loop counter. The value of abs(delta_i) could be saved in each loop and get compared to one from the previous iteration. If it is larger, i.e. the iteration is divergent, the program may exit with an error message. I am currently testing this in a 35s program.
Dieter
RE: 41C/CV root finders - Ángel Martin - 05-30-2015 02:04 PM
(05-30-2015 12:37 PM)Dieter Wrote: There is a similar method like the one suggested in post #13. The calculation of your f(i) and f'(i) could get implemented this way:
Code:
Let
a = expm1(-n * ln1p(i))
b = PMT * (1/i + p) - FV
Then
-f(i) = a * b - PV - FV
f'(i) = a * PMT/i^2 + b * n * (1+a)/(1+i)
Afterwards –f(i) is divided by f'(i) and this quotient is the delta_i that is then added (!) to the previous estimate. Now every iteration requires just one ln1+x and one e^x-1, the formulas are shorter and maybe even more accurate.
Thanks for the additional suggestions Dieter, indeed you're getting to the end of this and that's appreciated.
I tend to use more "linear" expressions in MCODE because they don't require any storage of intermediate values, like "a", and "b" above. Yes I know this saves time when reused in both formulas but the problem with 13-digit values is you run out of scratch registers very quickly.
Realize you'll need 4 standard registers to store a/b in 13-digit format, and that for convenience in the recall.store actions the stack is already used up with the TVM data. There's always the ALPHA registers of course, but I'm already using M,N, and O. I don't want to use any data register above R5 either.
Let me think about that, I like the approach and sure it's worth a try.
(05-30-2015 12:37 PM)Dieter Wrote: Regarding cases with no solutions and infinite looping: I think there is a better way than a simple loop counter. The value of abs(delta_i) could be saved in each loop and get compared to one from the previous iteration. If it is larger, i.e. the iteration is divergent, the program may exit with an error message. I am currently testing this in a 35s program.
Completely agree. Typically a diverging process won't recoup itself into a converging one, or will it?
'AM