2018年4月15日 星期日


Bisection Method To Solve Equation


Ching-Tang Tseng
Hamilton, New Zealand
16 April 2018

Catch up the personal computer developed performance. Using Bisection method to find real number roots of a given equation was getting into a very convenient condition. Originally, owing to the computer plotting feature is not so well, this numerical method in mathematics was not so welcome in public. For the present, some people still ignored its convenience and said that the absolute error in Bisection method is halved at each step so it converges linearly, which is comparatively slow.

Nowadays, personal computer running speed is fast enough. Speed is no longer to be a must-concerned problem. Contrarily, convenient feature let me tested out that other methods are usually lacked of absolute convergent ability but Bisection method has. Mathematics is a rigorous science. Divergent problem or trapped into a looping iteration problem makes many numerical techniques unable covered on full range of an application. At this point, Bisection method is the best.

Bisection method needs neither differential treatment nor some other complicated treatments. Its iteration is impossible getting into an infinite circulation condition. You are able to use this method to solve an equation over its whole definition domain. In this article, an integrate Bisection method is applied to solve a little bit complicated equation. Convenience is especially showing.

1. Introduction

There are many resources on many websites talked about Bisection method. I preferred to cite Bisection method material from our famous FORTH contributor J. V. Nobles website:


Section 13. Non-trivial programming example” in his book, he wrote the algorithm:

\*
If we know that the roots are bracketed between xa and xb, and that f(xa)*f(xb) < 0 (at least 1 root lies in the interval) we take the next guess to be xp = (xa+xb)/2 .
We then evaluate the function at xp: fp = f(xp).If fa*fp > 0 we set xa = xp, else we set xb = xp. We repeat until the ends of the interval containing the root are sufficiently close together.
\*

Such an explanation is simple enough and is sound enough, enough to apply to write code out directly.

By the way, an intention from my mind is to encourage everybody to read this J. V. Nobles book.

2. Example equation

y = ln |140 / (80-x) | - ( (60+x) / 53.0516477 )

28 years ago (1990), this equation had ever been used to test my FORTH infix function design. 10 years ago (2008), this equation had ever been used to test my ABC FORTH. Now (2018), I would like to use it to test my Lina64 ABC FORTH again. Original equation had no absolute value expression inside. This time, I added an absolute value function.

For the purpose of convenience, I used Lina64 to manipulate a Python plotting program to show the figure of this equation. This posting doesn’t mean that I had designed a Python in my Lina64 FORTH system. But Lina64 FORTH is really able to run Python in FORTH environment.

Python code:

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-100, 150, 1000)
y = np.log(np.absolute(140/(80-x)))-((60+x)/53.0516477)

plt.figure(figsize=(8,5))
plt.title("Fig.y=ln|140/(80-x)|-((60+x)/53.0516477) Designer:CT Tseng")

plt.xlabel("x")
plt.ylabel("y")

plt.plot(x,y,color="red",linewidth=4)

plt.ylim(-2,10)

plt.grid()
plt.axhline(y=0,color="green",linewidth=1)
plt.axvline(x=0,color="green",linewidth=1)

plt.show()


































Fig. 1

From Fig. 1, you are able to find out that there was three roots existed. The most left one is sited near 50. The middle one is above +50, the most right one is close to +100.
These values are the data I needed to use in my program.

3. Code

\ Bisection method to find real roots of equation in Lina64 ABC forth
\ Equation: y=ln|140/(80-x)|-((60+x)/53.0516477)
\ Guessing range: root1[-100, -50], root2[50, 80], root3[80, 100].

\ ===============================================
\ 1. Variables declare heap
variable Counter
integer Steps
10 reals X Y XA XB XP FA FB FP EPS TEMP

\ ===============================================
\ 2. Input heap
\ Equation is defined here.
: FX
{{ Y = LN ( ABS ( 140.0 e 0 / ( 80.0 e 0 - X ) ) ) - ( ( 60.0 e 0 + X ) / 53.0516477 e 0 ) }} ;

\ guessing range of roots is defined here.
: root1
{{ XA = -100.0 e 0 }} {{ XB = -50.0 e 0 }} ;

: root2
{{ XA = 50.0 e 0 }} {{ XB = 80.0 e 0 }} ;
: root3
{{ XA = 80.0 e 0 }} {{ XB = 100.0 e 0 }} ;
\ All parts above are changeable setting.
\ ===============================================

\ ===============================================
\ 3. Main code heap
\ All parts below are fixed.

: Fbisection
0 COUNTER !
BEGIN
{{ X = XA }} FX {{ FA = Y }}
{{ X = XB }} FX {{ FB = Y }}
{{ EPS = ABS ( XA - XB ) - 1.0 E -12 }}
EPS F0>
WHILE
{{ XP = ( XA + XB ) / 2.0 e 0 }}
{{ X = XP }} FX {{ FP = Y }}
{{ TEMP = FA * FP }}
TEMP F0>
IF {{ XA = XP }} ELSE {{ XB = XP }} THEN
1 COUNTER +!
REPEAT
CR ." Counts = " COUNTER @ .
{{ TEMP = ( XA + XB ) / 2.0 e 0 }}
CR ." Solution = " TEMP F. cr
;

: Bbisection
BASIC
10 let steps = 0

20 LET { X = XA }
30 RUN FX
40 LET { FA = Y }
50 LET { X = XB }
60 RUN FX
70 LET { FB = Y }
80 LET { EPS = ABS ( XA - XB ) - 1.0 E -12 }
90 IF { EPS <= 0.0 e 0 } THEN 210

100 LET { XP = ( XA + XB ) / 2.0 e 0 }
110 LET { X = XP }
120 RUN FX
130 LET { FP = Y }
140 LET { TEMP = FA * FP }
150 IF { TEMP > 0.0 e 0 } THEN 180
160 LET { XB = XP }
170 GOTO 190
180 LET { XA = XP }
190 LET Steps = Steps + 1

200 GOTO -20

210 PRINT " Steps = " ; steps
220 LET { TEMP = ( XA + XB ) / 2.0 e 0 }
230 run cr ." Solution = " TEMP f. cr
240 END
;
\ ===============================================

\ ===============================================
\ 4. Manipulation heap

: mainF
BASIC
10 run root1
20 run Fbisection
30 run root2
40 run Fbisection
50 run root3
60 run Fbisection
100 end
;

: mainB
root1 Bbisection
root2 Bbisection
root3 Bbisection ;

\ Russia: Moscow Radio Station signed voice
: Signal
MSO T/2 MSO T/2 HDO 3T/2 MSI T/2 MLA T/2 MSI T/2 HDO T/2 HRE T/2 HDO 1T MSO 2T ;
: signed
CTONE ADAGIO Signal 500 ms Signal
s" DeLaswedge! This, Is, Moscow, Radio! Station! " girlsay ;

: main
cr ." ===== FORTH style program outcomes ===== " cr
s" FORTH style program outcomes " girlsay
mainF
cr cr ." ===== BASIC style program outcomes ===== " cr
s" BASIC style program outcomes " boysay
mainB
cr ." ============ the end =================== " cr
s" The end " boysay
signed ;

4. Outcome

AMDX86 ciforth 5.2.1
fload bisection.f
FA OK
main

===== FORTH style program outcomes =====

Counts = 46
Solution = -60.00000000

Counts = 45
Solution = 67.291782187

Counts = 45
Solution = 88.517724806


===== BASIC style program outcomes =====

Steps = 46
Solution = -60.00000000

Steps = 45
Solution = 67.291782187

Steps = 45
Solution = 88.517724806

============ the end ===================
OK

5. Discussion

To plot an equation quickly out is very convenient now on my personal computer. This convenience makes me easy to set up two guessing values for each real number root. I did all operation under Lina64 system only. That means I am able to adjust all data in the program quickly and to check the figure easily.

Looking at Fig. 1, we know that there are real roots existed, and no more no less, total are three roots over there. Without such a figure we are unable to get such a conclusion so convenient.

When I was doing testing. I found that there was a singularity point x = 80. If the design of F/ operation was not smart enough to cover the divided by zero condition, adoption of guessing point took this point, system will be halted on divide by zero error. This problem could be checked out from Fig. 1 directly. In my system, F/ was able to deal with / 0 condition. If you did not asked system to print out an infinite or an -infinite message, program must kept going to go.

STEPS and COUNTS showed that less than 50 times iteration needed only. Iteration steps depend on two factors. The first factor is the tolerance setting, the smaller tolerance the more iteration steps, in my program, tolerance value EPS(epsilon) is setting to 1.0e-12. The second factor is the span range between lower bound and upper bound points. In my program, they are XA and XB. The wider span ranges the more iteration steps. Does the computing time a problem to me? No, it doesn’t matter to me absolutely, and computing time is really short enough.

In fact, I am able to make the span range narrower by checking the figure of equation. Zooming the curve and shrinking the x-axis scale range can reach it. That means convenience of plotting certainly can improve the Bisection method convergent speed.

The purpose of designing ABC FORTH is let me apt to do all mathematical calculation. This program has showed this point out. Program could be thought as a universal code. Equation definition in the program was mere one line and its coding format was the same as formal expression in mathematics. Special flavor could be felt out from reading these codes.

In the program, mainF is a BASIC style code. It calls FORTH style code. Correspond to this format, mainB is a FORTH style code. It calls BASIC style code.

My Lina64 has been accumulated up to 255 KB. This FORTH is a system: able to speak, able to sing, able to run many languages, Python and Perl included. A practiced Bisection method Perl program showed a feature that other programming language could be run by FORTH. Code is posting as following. Proof is showing in the video.

# Bisection method to find a root
#function defined
sub f
{
my $x = shift;
return ((60 + $x )/53.0516477)-log(140/(80 - $x));
}
#specified range and lower limit defined
my $xa = 60;
my $xb = 70;
my $epsilon = 1.0E-10;
#main program
while (abs($xa-$xb)>$epsilon)
{
my $x = $xa;
my $fa = &f($x);
my $xp = ($xa+$xb)/2;
my $x = $xp;
my $fp = &f($x);
if ($fa*$fp>0)
{
$xa = $xp;
}
else
{
$xb = $xp;
}
}
#output result
my $root = ($xa+$xb)/2;
print "root = $root\n";
#The end

Please appreciate the video demonstration at the end of this posting. Sometimes, I wrote my program is the same as to write an article.

6. Conclusion

In my personal opinion, Bisection method should be one subject of Must to learn list in numerical analysis. No matter what kind of the programming language student is learning, student ought to learn how to plot equations curve by program? Plotting and Bisection method combined together can solve out all real roots of an equation. Dont accept any criticism from an inexperienced talker, should listen to a good experience from a doer. Implement your personal code is the only one way to understand the real meaning of a Bisection method.

7. Video demonstration





8. My personal testing remarks

divided by 0 exception code in F/,
integrate Bisection method library code
approved.