HP50G Bisection methode

07242015, 12:01 PM
(This post was last modified: 07312015 03:51 PM by tigger.)
Post: #1




HP50G Bisection methode
I would like to program the Bisection Methode.
Is there anybody who can help me with it. I am an absolute beginner. I have the function y=x^55. I assume the result between 1 and 2: We first define the function F(X). Now save 2 in a variable called B and 1 in a variable called A. We have to find the midpoint between A and B, our new approximation, C. Can anbody help me to go on. I would like to program other easy stuff. 

07242015, 04:35 PM
(This post was last modified: 07242015 04:43 PM by Gilles.)
Post: #2




RE: HP50G Bisection methode
(07242015 12:01 PM)tigger Wrote: I would like to program the Bisection Methode.Hi It seems there is a problem with your example as there is no root for this function between 1 and 2 here is an example (not optimised), in RPL, aprox mode with a modified function (X^55) : Code: 'F(X)=X^55' DEFINE 

07292015, 11:31 AM
Post: #3




RE: HP50G Bisection methode
(07242015 12:01 PM)tigger Wrote: I would like to program the Bisection Methode. The root for y=x^5 (or for any y=x^n) is zero, not in the interval [1, 2]. You can start with an interval of, say, [1, 1]. Namir 

07292015, 11:37 AM
(This post was last modified: 07292015 11:39 AM by Namir.)
Post: #4




RE: HP50G Bisection methode
(07242015 04:35 PM)Gilles Wrote:(07242015 12:01 PM)tigger Wrote: I would like to program the Bisection Methode.Hi Your code assumes F(A)<0 and that's risky, because if F(A)>0 your implementation will not work. The traditional Bisection algorithm uses the test F(A)*F(C)>0 instead to determine if we set A=C (when that condition is true) and B=C otherwise. I assume the RPL code should be: Code: 'F(X)=X^55' DEFINE Namir 

10022020, 11:30 PM
(This post was last modified: 10032020 03:01 AM by acser.)
Post: #5




RE: HP50G Bisection methode
I was not able to get the code posted here earlier to work.
Checked https://www.math.ubc.ca/~pwalls/mathpyt...bisection/ for the algorithm of bisection. This worked for me: << > x << x 5 ^ 5  >> >> ‘F’ STO [Enter] 1 [Enter] 2 [Enter] << 0 > A B C << DO A B + 2. / DUP ‘C’ STO IF ‘F(A)*F(C)<0.’ THEN ‘B’ STO END IF ‘F(B)*F(C)<0.’ THEN ‘A’ STO END UNTIL ‘ABS(AB)<0.00000001’ END C >> >> [Enter] 

10032020, 02:03 AM
Post: #6




RE: HP50G Bisection methode
Hi, acser.
Welcome to the forum. Your code had a few bugs. Example, what happens if F(C) = 0 ? Also, we can reduce 4 functions calls per loop down to 1 Bonus: this avoided underflow issue, when F(A)*F(C) = 0.0, losing sign info. Code: def solve_bisect(f,a,b, eps=1e8, verbal=True): >>> f = lambda x: x**5  5 >>> solve_bisect(f,0,1) Traceback (most recent call last): File "<stdin>", line 1, in <module> File "<stdin>", line 6, in solve_bisect AssertionError: bad interval >>> solve_bisect(f,1,2) 1 2 1 1.5 1.25 1.5 1.375 1.5 1.375 1.4375 1.375 1.40625 1.375 1.390625 1.375 1.3828125 1.37890625 1.3828125 ... 

10042020, 03:04 PM
Post: #7




RE: HP50G Bisection methode
(10032020 02:03 AM)Albert Chan Wrote: >>> f = lambda x: x**5  5 For ε=1e15, solve_bisect() required log2(1/ε) ≈ 50 bisects (assuming it did not hit the root head on). I just learned Ridder's method. It has the advantage of bisection (guaranteed convergence), worstcase convergence of bisection. But, if function is well behaved, each loop approximately doubled convergence. (since it take 2 function calls to do 1 loop, convergence is √2 order) All code before getting (c, fc) is the same as solve_bisect(), then we add point (d, fd): Code: def solve_ridder(f,a,b, eps=1e8, verbal=True): >>> f = lambda x: x**5  5 >>> solve_ridder(f, 1,2, eps=1e15) 1 2 1.3789222674 1.5 1.37972953756 1.4394611337 1.37972966146 1.40959533563 1.37972966146 1.37972966146 1.37972966146 1.37972966146 1.3797296614612149 

10072020, 09:15 PM
(This post was last modified: 10072020 09:34 PM by acser.)
Post: #8




RE: HP50G Bisection methode
Thanks, Albert. I have not yet looked at the Ridders' method.
Here's the updated HP50g RPL code: Copy and paste the below code under the "" to the Emu48 stack then use the HP 50g program ASCIBIN.49 found at https://www.hpcalc.org/details/3648 to convert it to a proper HP 50g program on the stack. ‰ is the <= sign ([LS] [X]) and Š is the HP >= sign ([LS] [1/x])  %%HP: T(1)A(R)F(.); « IF OVER F OVER F > THEN SWAP END 0. 0. > R A B C FC « DO A B + 2. / 'C' STO 'F(C)' EVAL 'FC' STO IF FC 0. ‰ THEN C 'A' STO END IF FC 0. Š THEN C 'B' STO END UNTIL 'ABS(AB)< .0000001' END C » » 

10082020, 11:55 AM
Post: #9




RE: HP50G Bisection methode
(10072020 09:15 PM)acser Wrote: Copy and paste the below code under the "" to the Emu48 stack then use the HP 50g program ASCIBIN.49 found at https://www.hpcalc.org/details/3648 to convert it to a proper HP 50g program on the stack. Just as an aside to this interesting thread, the programs would be more usable if you used INOUT and put the code inside code tags. As it is, some comparisons in your program come out as odd special characters. 

10082020, 12:55 PM
Post: #10




RE: HP50G Bisection methode
Here you go:
'F(X)=X^55' DEFINE 1 [ENTER] 2 [ENTER] Use INOUT to get this code into your HP50g \<< IF OVER F OVER F > THEN SWAP END 0. 0. \> A B C FC \<< DO A B + 2. / 'C' STO 'F(C)' EVAL 'FC' STO IF FC 0. \<= THEN C 'A' STO END IF FC 0. \>= THEN C 'B' STO END UNTIL 'ABS(AB)<.0000001' END C \>> \>> 

« Next Oldest  Next Newest »

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