HP50G Bisection methode
07-24-2015, 12:01 PM (This post was last modified: 07-31-2015 03:51 PM by tigger.)
Post: #1
 tigger Junior Member Posts: 39 Joined: Jul 2015
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^5-5.
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.
07-24-2015, 04:35 PM (This post was last modified: 07-24-2015 04:43 PM by Gilles.)
Post: #2
 Gilles Member Posts: 171 Joined: Oct 2014
RE: HP50G Bisection methode
(07-24-2015 12:01 PM)tigger Wrote:  I would like to program the Bisection Methode.
Is there anybody who can help me with it.
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^5-5) :

Code:
'F(X)=X^5-5' DEFINE  1. 2.   «   IF OVER F OVER F > THEN SWAP END  0. -> A B C  «    DO    A B + 2. / DUP 'C' STO      IF 'F(C)<0' THEN 'A' STO ELSE 'B' STO END   UNTIL 'ABS(F(C))<0.01' END   C  » »
07-29-2015, 11:31 AM
Post: #3 Namir Senior Member Posts: 822 Joined: Dec 2013
RE: HP50G Bisection methode
(07-24-2015 12:01 PM)tigger Wrote:  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^5.
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.

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
07-29-2015, 11:37 AM (This post was last modified: 07-29-2015 11:39 AM by Namir.)
Post: #4 Namir Senior Member Posts: 822 Joined: Dec 2013
RE: HP50G Bisection methode
(07-24-2015 04:35 PM)Gilles Wrote:
(07-24-2015 12:01 PM)tigger Wrote:  I would like to program the Bisection Methode.
Is there anybody who can help me with it.
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^5-5) :

Code:
'F(X)=X^5-5' DEFINE  1. 2.   «   IF OVER F OVER F > THEN SWAP END  0. -> A B C  «    DO    A B + 2. / DUP 'C' STO      IF 'F(C)<0' THEN 'A' STO ELSE 'B' STO END   UNTIL 'ABS(F(C))<0.01' END   C  » »

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^5-5' DEFINE  1. 2.   «   IF OVER F OVER F > THEN SWAP END  0. -> A B C  «    DO    A B + 2. / DUP 'C' STO      IF 'F(A)*F(C)>0' THEN 'A' STO ELSE 'B' STO END   UNTIL 'ABS(F(C))<0.01' END   C  » »

Namir
10-02-2020, 11:30 PM (This post was last modified: 10-03-2020 03:01 AM by acser.)
Post: #5
 acser Junior Member Posts: 16 Joined: Sep 2020
RE: HP50G Bisection methode
I was not able to get the code posted here earlier to work.

Checked https://www.math.ubc.ca/~pwalls/math-pyt...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(A-B)<0.00000001’
END
C
>>
>>

[Enter]
10-03-2020, 02:03 AM
Post: #6
 Albert Chan Senior Member Posts: 1,791 Joined: Jul 2018
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=1e-8, verbal=True):     fa, fb = f(a), f(b)     if fa==0: return a     if fb==0: return b     if fa > fb: a, fa, b, fb = b, fb, a, fa     assert fa < 0 and fb > 0, "bad interval"     while abs(a-b) > eps:         if verbal: print a, b         c = (a+b)/2         fc = f(c)         if fc <= 0: a=c         if fc >= 0: b=c     return (a+b)/2

>>> 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

>>> 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
...
10-04-2020, 03:04 PM
Post: #7
 Albert Chan Senior Member Posts: 1,791 Joined: Jul 2018
RE: HP50G Bisection methode
(10-03-2020 02:03 AM)Albert Chan Wrote:  >>> f = lambda x: x**5 - 5
>>> 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
...

For ε=1e-15, 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), worst-case 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=1e-8, verbal=True):     fa, fb = f(a), f(b)     if fa==0: return a     if fb==0: return b     if fa > fb: a, fa, b, fb = b, fb, a, fa     assert fa < 0 and fb > 0, "bad interval"     while abs(a-b) > eps:         if verbal: print a, b         c = (a+b)/2         fc = f(c)         if fc == 0: return c         d = c - (c-a) * fc/(fc*fc-fa*fb)**0.5         fd = f(d)         if fd == 0: return d         if fd < 0:             a, fa = d, fd             if fc > 0: b, fb = c, fc         else:             b, fb = d, fd             if fc < 0: a, fa = c, fc     return (a+b)/2

>>> f = lambda x: x**5 - 5
>>> solve_ridder(f, 1,2, eps=1e-15)
1 2
1.3789222674 1.5
1.37972953756 1.4394611337
1.37972966146 1.40959533563
1.37972966146 1.37972966146
1.37972966146 1.37972966146
1.3797296614612149
10-07-2020, 09:15 PM (This post was last modified: 10-07-2020 09:34 PM by acser.)
Post: #8
 acser Junior Member Posts: 16 Joined: Sep 2020
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(A-B)<
.0000001'
END C
»
»
10-08-2020, 11:55 AM
Post: #9 John Keith Senior Member Posts: 751 Joined: Dec 2013
RE: HP50G Bisection methode
(10-07-2020 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.

‰ 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(A-B)<
.0000001'
END C
»
»

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.
10-08-2020, 12:55 PM
Post: #10
 acser Junior Member Posts: 16 Joined: Sep 2020
RE: HP50G Bisection methode
Here you go:

'F(X)=X^5-5' 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(A-B)<.0000001'
END C
\>>
\>>
 « Next Oldest | Next Newest »

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