Wednesday 12 October 2011

Root finding using the Bisection Method - comments on the logic

In my previous blog entry I described one way in which to implement the bisection method. Most programmers would probably be aghast at my liberal use of the Goto command, but in this situation I think it is required.

There are necessarily two calculations of the function f(x) in different places of the logic; f(a) at the beginning of the program and f(c) during the iterations. As I only want to have the function in one place in the program (to make it easier to edit) I needed to have a subroutine, but the fx-50F doesn't have this functionality. To make matters worse, the control structures such as While/WhileEnd do not allow you to jump in and out of them using a Goto command and they also cannot be nested. So in the end I had to go back to basics and use Goto commands with a structure similar to a subroutine and a While/WhileEnd loop.

The 'subroutine' can be found at the end of the logic after the statement Display X (this is the actual end-point of the program). I have included another memory variable X which counts the number of iterations that were needed to obtain the root to the required accuracy. Conveniently, X has been used to identify where in the program the subroutine has been called from.

The function that appears in this subroutine is f(x)=x-cos(x) and I will say more about this in a later blog.

Wednesday 5 October 2011

Root finding using the Bisection Method - the logic

We want to come up with an algorithm for finding the root of a function f(x) at x=α (where f(α)=0) using the bisetion method. We have 6 key pieces of information which we need to store in the calculator's memory as the iteration process described takes place. For the interval a<x<b being bisected, these are the value of a, b, δ, f(a), c=(a+b)/2, and f(c). Let the memory variables which store these values be A, B, D, M, C and Y, respectively. δ is the accuracy to which we wish to calculate α.

On the first iteration we key in the values a and b which span the root and store them in A and B. We also key in the value of δ and save it in D. We calculate c=(a+b)/2, f(a) and f(c) and store them in C, M and Y, respectively. Our estimate of the value of α in this first iteration is c and the error in the estimate is (b-a)/2. So if (b-a)/2 is less than our tolerance δ or if f(c)=0 then we are finished and we can display the contents of C as our estimate of α.

If, on the other hand, (b-a)/2 is greater than our tolerance δ and f(c) is not zero we can determine whether the root lies in the range a<x<c or c<x<b. If the product f(a)f(c)>0 then the root lies in the interval c<x<b otherwise it lies in the interval a<x<c. So if the product MY>0 we can store C in A and Y in M otherwise we store C in B. This gives us our new end points a and b (held in A and B) to process in the next iteration and, possibly, a new value of f(a) in M.

On the second and subsequent iterations we recalculate c=(a+b)/2 and f(c) using the values in A and B and store them in C and Y, respectively. We can then repeat the processing that we did in the first iteration.

The processing will look something like this:-

Input A
Store A in C
Store 0 in X
Goto Label 1
Label 2
Store Y in M
Input B
Input D
Label 4
Store X+1 in X
Store 0.5(A+B) in C
If 0.5(B-A) is less than D then Goto Label 0
Goto Label 1
Label 3
If Y=0 then Goto Label 0
If MY is greater than 0 then
   Store C in A
   Store Y in M
Else
   Store C in B
IfEnd
Goto Label 4
Label 0
Display C
Display X
Label 1
Store C-Cos(C) in Y
If X=0 then Goto Label 2
If X is not equal to 0 Goto Label 3