Post Reply 
HP50G Bisection methode
07-24-2015, 12:01 PM (This post was last modified: 07-31-2015 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^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.
Find all posts by this user
Quote this message in a reply
07-24-2015, 04:35 PM (This post was last modified: 07-24-2015 04:43 PM by Gilles.)
Post: #2
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
 »
»
Find all posts by this user
Quote this message in a reply
07-29-2015, 11:31 AM
Post: #3
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
Find all posts by this user
Quote this message in a reply
07-29-2015, 11:37 AM (This post was last modified: 07-29-2015 11:39 AM by Namir.)
Post: #4
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
Find all posts by this user
Quote this message in a reply
10-02-2020, 11:30 PM (This post was last modified: 10-03-2020 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/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]
Find all posts by this user
Quote this message in a reply
10-03-2020, 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=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
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
...
Find all posts by this user
Quote this message in a reply
10-04-2020, 03:04 PM
Post: #7
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
Find all posts by this user
Quote this message in a reply
10-07-2020, 09:15 PM (This post was last modified: 10-07-2020 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(A-B)<
.0000001'
END C
»
»
Find all posts by this user
Quote this message in a reply
10-08-2020, 11:55 AM
Post: #9
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.
Find all posts by this user
Quote this message in a reply
10-08-2020, 12:55 PM
Post: #10
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
\>>
\>>
Find all posts by this user
Quote this message in a reply
Post Reply 




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