A somewhat different Newton solver (HP35s)

08072016, 08:47 PM
(This post was last modified: 08072016 08:48 PM by Dieter.)
Post: #1




A somewhat different Newton solver (HP35s)
I think most of us will know the usual Newton method for finding the roots of a function. This requires the evaluation of the function f(x) as well as its derivative f'(x), which usually means two function calls per iteration step.
For a project where the derivative cannot be calculated fast and easily I now tried the following method that is based on a method already discussed earlier, e.g. here and there. The idea here is the simulaneous calculation of f(x) and f'(x) via complex math. On the 35s this can be done quite easily – at least within the supported function set. Here is a sample program to illustrate the idea. Code: S001 LBL S Please note: this is just to show the idea, not a fully working program. ;) Steps S009...S014 calculate the quotient of the real (approx. f(x)) and imaginary parf (approx. f'(x)). The sine and cosine should not be replaced by 1/tangent since f(x)=0 will yield a real part of ±90° where the tangent is not defined (a cotangent function would be nice here...). Example: Code your function f(x) at LBL F, assuming the argument x in the Xregister. The sample function is x³–x²–x+0,5 = 0. After XEQ S enter an initial guess for x. Code: [XEQ] S [ENTER] X? So this seems to work. What about limitations of this method? As far as I get it both f(x) as well as f'(x) are only approximate. I wonder how accurate especially f(x) is evaluated using this method. Maybe the math experts here can say more about this. ;) Dieter 

09232018, 05:28 PM
(This post was last modified: 09242018 01:59 AM by Thomas Klemm.)
Post: #2




RE: A somewhat different Newton solver (HP35s)
(08072016 08:47 PM)Dieter Wrote: What about limitations of this method? As far as I get it both f(x) as well as f'(x) are only approximate. I wonder how accurate especially f(x) is evaluated using this method. Maybe the math experts here can say more about this. ;) Using the Taylor series we can separate the real and imaginary part of \(f(x+ih)\): \(\begin{align*} \Re[f(x+ih)]&=f(x)\frac{f''(x)}{2!}h^2+O(h^4)\\ \Im[f(x+ih)]&=h(f'(x)\frac{f'''(x)}{3!}h^2+O(h^4)) \end{align*}\) Thus by choosing a small \(h\) we can make the difference as small as we want. You can verify this by using \(h=10^{50}\) and the \(\sin\) function: sin(2 + 1e50 i) = 9.09297426826e1  4.16146836547e51 i sin(2) = 9.09297426826e1 cos(2) = 4.16146836547e1 Both values are exact. But for a calculator with 12 significant digits something smaller like \(10^{10}\) might be enough. It really depends on the value of the 2nd and 3rd derivative at the root. We could go further and try to use \(h=10^{99}\). But then we may end up with 0. Thus we have to choose \(h\) small enough to avoid errors but not too small. Kind regards Thomas 

09252018, 11:46 AM
(This post was last modified: 09252018 01:54 PM by Ángel Martin.)
Post: #3




RE: A somewhat different Newton solver (HP35s)
Nicely done, it piqued my curiosity so I went ahead and jotted down the following 41Z version of the same concept. About 30 steps in this proof of concept; which expects the function name in ALPHA, h in Y and the xo guess (a real value) in X.
Code: 01 LBL "ZNWT" The execution shows the successive values of the root if user flag 10 is set (a la PPCROM). Convergence criteria always uses 9 decimal places  irrespective of the display FIX settings. The function is to be programmed using 41Z functions (definitely a superset of those in the 35S), as a complex equation. To use Dieter's same example: Code: 01 LBL "Z1" Results: ALPHA, "Z1", ALPHA 0.1, ENTER^, 2, XEQ "ZNWT" > 1.451605963 0.1, ENTER^, 0, XEQ "ZNWT" > 0.403031717 0.1, ENTER^, 2, XEQ "ZNWT"> 0.854637680 Supposedly more accurate results would be obtained with smaller values of h... Cheers, ÁM 

09262018, 07:56 AM
Post: #4




RE: A somewhat different Newton solver (HP35s)
You don't have to execute the function twice since both, value and derivative are returned:
Code: 01 LBL "ZNWT" ; h x Example: 2 ENTER 1E9 XEQ "ZNWT" 2.000000000 1.642857143 1.487473705 1.453261463 1.451609751 1.451605963 (09252018 11:46 AM)Ángel Martin Wrote: The execution shows the successive values of the root if user flag 10 is set (a la PPCROM). Feel free to change that to your needs: Add the check of user flag 10 and remove RND. Quote:Supposedly more accurate results would be obtained with smaller values of h... Usually 1e9 will be small enough. But you might also use 1e50. Kind regards Thomas 

09262018, 09:05 AM
(This post was last modified: 09262018 11:21 AM by Ángel Martin.)
Post: #5




RE: A somewhat different Newton solver (HP35s)
(09262018 07:56 AM)Thomas Klemm Wrote: You don't have to execute the function twice since both, value and derivative are returned: Brilliant Thomas, many thanks for snapping me out of my stupor to see what was before my nose! And so glad to see there's at least one more 41Zliterate person in the world ;) BTW nice move choosing R02 for the function name instead of R00  this allows using the 41Z functions on ZR00, which uses "0" as default parameter so it's not required in the program. Really impressive  although you must be using an older revision of the module, because with the latest one (the 'Deluxe" edition) the real part is stored in the lower register of the pair, i.e. ZSTO 00 stores Re(z) in R00 and Im(z) in R01 So rewriting the code for the 41ZDeluxe conventions (transposing R00 and R01), and starting with h in Y, x guess in X: Code: 01 LBL "ZNWT" ; x h 17 steps, nothing can beat that ;) Thanks much again. Best, ÁM PS. As it turns out removing RND is not a good idea because some times the oscillation in the 10th. decimal place causes the test X#0? never to occur. This is a wellknown behavior of the Newton method, so nothing new here... 

09262018, 12:49 PM
Post: #6




RE: A somewhat different Newton solver (HP35s)
(09262018 09:05 AM)Ángel Martin Wrote: And so glad to see there's at least one more 41Zliterate person in the world ;) You made me read your excellent documentation. And who in this forum doesn't love to read manuals? Quote:BTW nice move choosing R02 for the function name instead of R00  this allows using the 41Z functions on ZR00, which uses "0" as default parameter so it's not required in the program. I've tried that but somehow it didn't work. But I must admit I didn't invest too much into it. Quote:Really impressive  although you must be using an older revision of the module, because with the latest one (the 'Deluxe" edition) the real part is stored in the lower register of the pair, i.e. ZSTO 00 stores Re(z) in R00 and Im(z) in R01 I was using a somewhat buggy/unstable emulator my41CX on my iPhone since I didn't have access to an emulator on a computer. When the calculator is turned off it displays for a moment: 41ZREV: 9Z Cheers Thomas 

09272018, 12:06 PM
Post: #7




RE: A somewhat different Newton solver (HP35s)
(09262018 09:05 AM)Ángel Martin Wrote: I don't think it's a good idea to round the calculated ∆x to display precision and then adjust x by this lowerprecision value. So you should first adjust x (ST– 00) and then RND the delta and do the X≠0 test. Simply swap line 12 and 13. Dieter 

09272018, 12:38 PM
(This post was last modified: 09272018 12:40 PM by Albert Chan.)
Post: #8




RE: A somewhat different Newton solver (HP35s)
(08072016 08:47 PM)Dieter Wrote: I think most of us will know the usual Newton method for finding the roots of a function. This requires the evaluation of the function f(x) as well as its derivative f'(x), which usually means two function calls per iteration step. Does the project reach its goal ? I know complex math can get better derivative, but is it really fast ? I would guess complex math for f(x) and f(x + h*I) are much slower than real f(x) and f(x+h). Just look at effort needed to multiply of 2 complex numbers ... 

09272018, 12:39 PM
(This post was last modified: 09272018 12:44 PM by Ángel Martin.)
Post: #9




RE: A somewhat different Newton solver (HP35s)
(09272018 12:06 PM)Dieter Wrote: I don't think it's a good idea to round the calculated ∆x to display precision and then adjust x by this lowerprecision value. So you should first adjust x (ST– 00) and then RND the delta and do the X≠0 test. Simply swap line 12 and 13. Very true, thanks. (09262018 12:49 PM)Thomas Klemm Wrote: I was using a somewhat buggy/unstable emulator my41CX on my iPhone since I didn't have access to an emulator on a computer. Yes, that's a very old revision, probably not even using the Library#4... 

09272018, 12:41 PM
(This post was last modified: 09272018 12:45 PM by Ángel Martin.)
Post: #10




RE: A somewhat different Newton solver (HP35s)
(09272018 12:38 PM)Albert Chan Wrote: Does the project reach its goal ? I guess it all depends on the platform. On the 41Z it definitely meets the goal  with a single execution of the function to calculate both the function *and* its derivative. I'd think that's also the case for the 42S and even the 35S since all complex functions are also MCODE. 

09272018, 01:07 PM
(This post was last modified: 09272018 01:08 PM by Dieter.)
Post: #11




RE: A somewhat different Newton solver (HP35s)
(09272018 12:38 PM)Albert Chan Wrote: Does the project reach its goal ? A complex multiply is just four real multiplications, and two additions/subtractions. Thats "nothing". While I can't say much about this method in high level languages with complex number support by the compiler, I bet that on almost any calculator this method will be faster than the classic (approximative) approach with f(x) and f(x+h): First, there is just one single function call per iteration step. The complex result holds both f(x) and f'(x). The classic method requires two calls per iteration. This alone speeds up the calculation, even if the complex functions should be somewhat slower. More important, I don't know any programmable pocket calculator where complex functions are returned significantly slower to the user than their real counterparts. Sure, the internal calculations are more elaborate, but all this is "low level" code. The execution time for such functions is completely negligible compared to the running speed of the calculator program, which does merely a dozen or maybe a few hundred user instructions per second. The CPU that handles the complex math "under the hood" may do thousands of low level instructions per second to calculate a function that is returned almost immediately. So all in all the complex evaluation is not significantly slower than real math, and at the same time only half of the function calls is required. Yes, that's faster. ;) Dieter 

09272018, 10:04 PM
Post: #12




RE: A somewhat different Newton solver (HP35s)
.
Hi, Dieter: (09272018 01:07 PM)Dieter Wrote: A complex multiply is just four real multiplications, and two additions/subtractions. It can be done with just three real multiplications (instead of 4) and five additions/subtractions (instead of two), i.e. trading one multiplication for three add/subs. In user code (say RPN) this would probably fail to be advantageous because the time to execute a "*" is about the same as the time to execute a "+" or "", which in both cases is dominated by overhead not related to the actual math operation. However, if doing the multiplications, additions and subtractions as part of a microcoded keyword or function (such as the ones Ángel implements all the time) then the individual overheads are much less and it might be the case that calling (or implementing) a microcode (or assembly language) multiplication might be more expensive than doing three extra subs/adds, and thus the 3m+5as way would be faster than the 4m+2as one. Ángel probably has tested this possibility and if so he knows what's better (or if not he might try it) and will perhaps share his experience with us. Regards. V. . All My Articles & other Materials here: Valentin Albillo's HP Collection 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)