Der Zwischenwertsatz besagt: Jede stetige Funktion $f$ nimmt im Intervall $[a,b]$ jeden Wert von $f(a)$ und $f(b)$ mindestens einmal an. Bei einem Vorzeichenwechsel im Intervall $[a,b]$ muss es also mindestens eine Nullstelle geben.

Nicht das schnellste, aber das robustete Verfahren ist die Intervallhalbierung. Untersucht wird das Vorzeichen in der Mitte $m$ der Intervalls. Je nachdem, ob der Vorzeichenwechsel in der Hälfte $[a,m]$ oder in der Hälfte $[m,b]$ stattfindet, wird dort weiter gesucht, bis die Genauigkeitsanforderungen erfüllt sind oder eine genügende Anzahl von Halbierungsschritten unternommen wurde. Am Ende wird ein (engeres) Intervall geliefert, das die Nullstelle sicher einschließt:

bisektion.cpp

#include <cassert>
#include <cmath>
#include <iostream>
#include <utility>
 
template <typename Domain>
auto midpoint(Domain a, Domain b) { return (a + b) / 2; }
 
template <typename Function, typename Domain>
std::pair<Domain, Domain> 
root_bisection(Function f, Domain a, Domain b, int steps = 50)
{
  assert(f(a) * f(b) < 0);
  auto ya = f(a);
 
  while (steps-- > 0)
  {
    auto m = midpoint(a, b);
    auto ym = f(m);
    if (ya * ym < 0) { b = m; }
    else { a = m; ya = ym; }
  }
  return {a, b};
}
 
int main()
{
  auto f = [](double x) { return std::sin(x); };
  auto x = root_bisection(f, 3.0, 4.0);
 
  auto a = x.first, b = x.second;
  std::cout << a << " " << b <<  " : " << f(a) << " " << f(b) << '\n';
  return 0;
}

Die Boost.Math-Bibliothek enthält ähnliche Routinen, setzt aber reellwertige Funktionen mit reellem Definitionsbereich voraus, während hier auch andere Grundbereiche (Domain) erlaubt sind. Template-Varianten mit vorgegebenem Abstand distance(Domain a, Domain b) sind leicht abzuleiten und durch Spezialisierung der Funktionen midpoint(a,b) bzw. distance(a,b) erweiterbar:

bisektion.h

#include <cassert>
#include <cmath>
#include <utility>
 
template <typename Domain>
auto midpoint(Domain a, Domain b) { return (a + b) / 2; }
 
template <typename Domain>
auto distance(Domain a, Domain b) { using std::abs; return abs(a - b); }
 
template <typename Function, typename Domain, typename Distance>
std::pair<Domain, Domain> 
root_bisection(Function f, Domain a, Domain b, Distance epsilon, int steps = 50)
{
  assert(f(a) * f(b) < 0);
  auto ya = f(a);
 
  while (distance(a, b) > epsilon && steps-- > 0)
  {
    auto m = midpoint(a, b);
    auto ym = f(m);
    if (ya * ym < 0) { b = m; }
    else { a = m; ya = ym; }
  }
  return {a, b};
}