Threaded Mode | Linear Mode

tigger · (This post was last modified: 07-31-2015 03:51 PM by tigger.)

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.

Gilles · (This post was last modified: 07-24-2015 04:43 PM by Gilles.)

(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

 »

»

Namir · 07-29-2015, 11:31 AM

(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

Namir · (This post was last modified: 07-29-2015 11:39 AM by Namir.)

(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

acser · (This post was last modified: 10-03-2020 03:01 AM by acser.)

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]

Albert Chan · 10-03-2020, 02:03 AM

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

Albert Chan · 10-04-2020, 03:04 PM

(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

acser · (This post was last modified: 10-07-2020 09:34 PM by acser.)

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

John Keith · 10-08-2020, 11:55 AM

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

acser · 10-08-2020, 12:55 PM

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