Monday, March 28, 2011

Bisection Method

Pseudo-code

                            Here is a representation of the bisection method in Pseudo-code. The variables left and right correspond to a and b above. The initial left and right must be chosen so that f(left) and f(right) are of opposite sign (they 'bracket' a root). The variable epsilon specifies how precise the result will be.
'Bisection Method
 
 'Start loop
 Do While (abs(right - left) > 2*epsilon)
 
   'Calculate midpoint of domain
   midpoint = (right + left) / 2
 
   'Find f(midpoint)
   If ((f(left) * f(midpoint)) < 0) Then
     'Throw away right half
     right = midpoint
   ElseIf ((f(right) * f(midpoint)) < 0)
     'Throw away left half
     left = midpoint
   Else
     'Our midpoint is exactly on the root
     Exit Do
   End If
 Loop
 Return midpoint

[edit] Practical considerations

For robust usage in practice, some care is required due to the properties of floating-point arithmetic.
  1. The expression f(left) * f(midpoint) is very likely to underflow to 0 since both arguments are approaching a zero of f. To avoid this possibility, the two function values should be tested separately rather than multiplied.
  2. If epsilon is too small, the value of abs(right - left) might never become as small as 2*epsilon, as left and right can get stuck at adjacent non-equal floating-point values. This possibility can be avoided by disallowing epsilon to be too small (depending on the precision of the arithmetic) or by adding extra tests to detect the stuck condition. Another method relies on the assumption that (right+left)/2 exactly equals left or right if left and right are adjacent non-equal floating-point values. This is true for most hardware, including IEEE arithmetic in the absence of overflow, but can be violated if intermediate expressions are calculated to greater precision than stored variables (depending on computer optimization).
Code which takes both these problems into account (and relies on the uncertain assumption) is as follows. It attempts to find a zero as accurately as the machine precision allows.
// Assumption: One of f(a) and f(b) is ≥ 0 and the other is ≤ 0
  if f(a) <= 0 then
      lo := a; hi := b
  else
      lo := b; hi := a
  endif
 
  mid := lo + (hi-lo)/2
  while (mid ≠Pseudo-codeHere is a representation of the bisection method in Pseudo-code. The variables left and right correspond to a and b above. The initial left and right must be chosen so that f(left) and f(right) are of opposite sign (they 'bracket' a root). The variable epsilon specifies how precise the result will be.
'Bisection Method
 
 'Start loop
 Do While (abs(right - left) > 2*epsilon)
 
   'Calculate midpoint of domain
   midpoint = (right + left) / 2
 
   'Find f(midpoint)
   If ((f(left) * f(midpoint)) < 0) Then
     'Throw away right half
     right = midpoint
   ElseIf ((f(right) * f(midpoint)) < 0)
     'Throw away left half
     left = midpoint
   Else
     'Our midpoint is exactly on the root
     Exit Do
   End If
 Loop
 Return midpoint

[edit] Practical considerations

For robust usage in practice, some care is required due to the properties of floating-point arithmetic.
  1. The expression f(left) * f(midpoint) is very likely to underflow to 0 since both arguments are approaching a zero of f. To avoid this possibility, the two function values should be tested separately rather than multiplied.
  2. If epsilon is too small, the value of abs(right - left) might never become as small as 2*epsilon, as left and right can get stuck at adjacent non-equal floating-point values. This possibility can be avoided by disallowing epsilon to be too small (depending on the precision of the arithmetic) or by adding extra tests to detect the stuck condition. Another method relies on the assumption that (right+left)/2 exactly equals left or right if left and right are adjacent non-equal floating-point values. This is true for most hardware, including IEEE arithmetic in the absence of overflow, but can be violated if intermediate expressions are calculated to greater precision than stored variables (depending on computer optimization).
Code which takes both these problems into account (and relies on the uncertain assumption) is as follows. It attempts to find a zero as accurately as the machine precision allows.
// Assumption: One of f(a) and f(b) is ≥ 0 and the other is ≤ 0
  if f(a) <= 0 then
      lo := a; hi := b
  else
      lo := b; hi := a
  endif
 
  mid := lo + (hi-lo)/2
  while (mid ≠ lo) and (mid ≠ hi) do
     if f(mid) ≤ 0 then
        lo := mid
     else
        hi := mid
     endif
     mid := lo + (hi-lo)/2
  endwhile
 
  return mid
lo) and (mid ≠ hi) do if f(mid) ≤ 0 then lo := mid else hi := mid endif mid := lo + (hi-lo)/2 endwhile return mid