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. Noble’s
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. Noble’s
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 equation’s
curve by program? Plotting and Bisection method combined together
can solve
out all real roots of an equation. Don’t
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.
沒有留言:
張貼留言