FORTRAN Help: Bisection Method & Roots of Functions

[mattmac.nuke]

Homework Statement

The purpose of this program is to calculate the approximate roots of the Sine function on given intervals. The intervals are input by the user, and then the do loop continues until the condition (m becomes very close to 0 or equals 0) is met.

The Attempt at a Solution

program bisec
IMPLICIT NONE
REAL :: a, b, m, f_xa, f_xb, f_xm

WRITE (*,*) 'Please enter the interval [A,B]:'
READ (*,*) a,b

DO !WHILE (ABS(m) > 1E-7)
m = (a + b)/2.
f_xa = SIN(a)
f_xb = SIN(b)
f_xm = SIN(m)

IF (ABS(m) < 1E-5) THEN
EXIT
END IF

IF (f_xa*f_xm > 0) THEN
a= m
ELSE IF (f_xa*f_xm < 0) THEN
b= m
ELSE IF (f_xa*f_xm == 0) THEN
EXIT
END IF

END DO
WRITE (*,*) 'Solution is:',m

end program bisec

This is what I have so far, I've changed around my conditional statement to see if it would help, but it did not. The solutions it gives me are either 0 (on an interval that does not include 0) or radically large numbers, or it runs indefinitely. I know I am doing something wrong within the do loop. I've tried to follow the math correctly and translate it into code, but this is still a challenge for me as I am still in the early learning stages of programming. Thank you!

 

[gneill]

What is the value of m the first time the DO condition is tested (that is, what is its value just prior to entering the DO loop for the first time?

Surely if you're looking for a root you want the value of the function f(m) to approach zero, not the value of m; m is simply the midpoint x-value, which takes on values in your search interval [a,b].

I note that you don't do any preliminary testing to see that the interval actually contains a zero. This may cause you grief (like an infinite loop for some pairs [a,b]).

 

[mattmac.nuke]

I made some modifications last night, after I posted actually. You're right, I'm not sure why on Earth I was test m instead of f_xm for the condition, but I've fixed that. I wound up initializing m inside the do loop, so it's initial value is the result of the operations on A and B, and it updates through each iteration. I thought about testing end points too, the program seems to work fine for positive and negative values where one of the limits (A or B) is actually the solution. The problem I seem to be having with negative intervals though is that when I punch in evaluate from [-8,-7] it will display -7 as a solution, however if I enter it as [-7,-8] it will display -8 as the correct solution. When evaluating for positive values however, no matter what the order of the numbers, it displays the correct solution. Why does it not do this for negative values as well?
Here is my updated code:

program bisec
IMPLICIT NONE
REAL(10) :: a, b, m, f_xa, f_xb, f_xm, pi
INTEGER :: i, n
WRITE (*,*) 'LAB 5- Bisection Method'
WRITE (*,*) '-----------------------'
WRITE (*,*) ' '
WRITE (*,*) 'Please enter the interval [A,B]:'
READ (*,*) a,b
pi = 3.1415926

DO i = 0, 1000
m = (a + b)/2.
f_xa = SIN(pi*a/2.)
f_xb = SIN(pi*b/2.)
f_xm = SIN(pi*m/2.)

IF (f_xa*f_xm > 0) THEN
a= m
ELSE IF (f_xa*f_xm < 0) THEN
b= m
ELSE IF (f_xa*f_xm == 0) THEN
EXIT
END IF

IF (ABS(f_xm) < 0.000000000000001) EXIT

END DO
WRITE (*,'(A, F8.3)') 'The solution on this interval is: x =',m
end program bisec


I really appreciate the assistance, Thank you!

 

[gneill]

Why don't you play computer and step, by hand and hand calculator, through one or two loops for one of the 'problematical' cases? Or work from a sketch of sin(x), since you really only need to verify the conditional logic for the choice of which end-point to move to 'm'.