Metoda Brenta

W analizie numerycznej metoda Brenta jest metodą znajdowania pierwiastków, która łączy metodę bisekcji, metodę siecznych oraz odwrotną interpolację kwadratową. Ma niezawodność bisekcji, ale może być tak szybka jak te mniej wiarygodne metody. Algorytm próbuje użyć potencjalnie szybko zbieżnych metod jak metodę siecznych czy, jeśli to możliwe, odwrotną interpolację kwadratową, ale wraca do bardziej solidnej metody bisekcji, jeśli to konieczne. Metodę tę wymyślił Richard Brent, opierając się na wcześniejszym algorytmie Theodorusa Dekkera. W związku z tym metoda jest znana pod nazwą metody Brenta-Dekkera.

Opis

Metoda Dekkera sprawuje się dobrze, jeśli funkcja $f$ stosunkowo dobrze się zachowuje. Jednak istnieją okoliczności, w których każda iteracja wykorzystuje metodę siecznych, ale iteracje $b_{k}$ zbiegają bardzo powoli (w szczególności $|b_{k}-b_{k-1}|$ może być dowolnie małe). W tym wypadku metoda Dekkera wymaga znacznie więcej iteracji niż metoda bisekcji. Brent w 1973 zaproponował małą modyfikację, aby uniknąć tego problemu. Wstawił dodatkowy test, który musi być spełniony zanim w następnej iteracji zostanie zaakceptowany rezultat metody siecznych. Muszą być spełnione jednocześnie dwie nierówności: biorąc pod uwagę właściwą tolerancję numeryczną, jeśli poprzedni krok użył metody bisekcji, musi zostać spełniona nierówność

$|\delta |<|b_{k}-b_{k-1}|,$

aby przeprowadzić interpolację. W przeciwnym przypadku wykonywana jest metoda bisekcji i jej rezultat jest użyty w następnej iteracji. Jeśli poprzedni krok wykonał interpolację, wtedy zamiast niej użyta jest nierówność

$|\delta |<|b_{k-1}-b_{k-2}|,$

do wykonania następnej akcji (do wyboru): interpolacji (kiedy nierówność jest spełniona) lub bisekcji (kiedy nierówność nie jest spełniona). Również, jeśli w poprzednim kroku została użyta metoda bisekcji, musi być spełniona nierówność

$|s-b_{k}|<{\tfrac {1}{2}}|b_{k}-b_{k-1}|.$

W przeciwnym razie wykonywana jest metoda bisekcji i jej rezultat jest użyty w następnej iteracji. Jeśli natomiast w poprzednim kroku użyto interpolacji, zamiast tego użyta jest nierówność:

$|s-b_{k}|<{\tfrac {1}{2}}|b_{k-1}-b_{k-2}|.$

Modyfikacja powoduje pewność, że na k-tej iteracji krok bisekcji zostanie wykonany w najwyżej $2\log _{2}(|b_{k-1}-b_{k-2}|/\delta )$ dodatkowych iteracjach, ponieważ powyższe warunki zmniejszają o połowę rozmiar kolejnego kroku interpolacji przy każdych dwóch iteracjach i po najwyżej $2\log _{2}(|b_{k-1}-b_{k-2}|/\delta )$ iteracjach wielkość kroku będzie mniejsza niż $\delta ,$ które wywołuje krok bisekcji. Brent dowiódł, że jego metoda wymaga co najwyżej N² iteracji, gdzie N oznacza ilość iteracji metody bisekcji. Jeśli funkcja $f$ dobrze się zachowuje, wtedy metoda Brenta postępuje za pomocą odwrotnie kwadratowej bądź liniowej interpolacji, co będzie powodować zbieżność ponadliniową. Co więcej, metoda Brenta używa odwrotnej interpolacji kwadratowej zamiast liniowej interpolacji (która jest użyta w metodzie siecznych). Jeśli $f(b_{k}),$ $f(a_{k})$ i $f(b_{k-1})$ będą różne, to wydajność zostanie nieznacznie zwiększona. W konsekwencji następuje zmiana warunku do zaakceptowania s (wartość proponowana zarówno przez liniową interpolację, jak i odwrotną interpolację kwadratową): s powinno leżeć pomiędzy $(3a_{k}+b_{k})/4$ i $b_{k}.$

Algorytm

Według oryginalnego programu Brenta w Algolu z poprawkami.

Brent(a, b, t) //t = zadana tolerancja, która będzie powiększana o błąd maszynowy
  {
      wylicz fa = f(a);
      wylicz fb = f(b);
      c = a; fc = fa;
      d = e = b - a;
      while (true)
      {
          if (|fc| < |fb|)
          {
              SWAP(b,c);
              SWAP(fb,fc);
          }
          tol = 2 * delta(b) + t;
          m = 0.5*(c - b);
          if (|m| > tol && fb != 0)
          {
              if (|e| < tol || |fa| <= |fb|)
                  d = e = m; //bisekcja
              else
              {
                  s = fb / fa;
                  if (a == c)
                  { // metoda siecznych
                      p = 2 * m * s;
                      q = 1 - s;
                  }
                  else
                  { // odwrotna interpolacja kwadratowa
                      q = fa / fc;
                      r = fb / fc;
                      p = s*(2 * m*q*(q - r) - (b - a)*(r - 1));
                      q = (q - 1)*(r - 1)*(s - 1);
                  }
                  if (p > 0) q = -q; else p = -p;
                  s = e; e = d;
                  if (2 * p < 3 * m*q - |tol*q| && p < |0.5*s*q|)
                      d = p / q;
                  else
                      d = e = m; //bisekcja
              }
              a = b; fa = fb;
              if (|d| > tol)
                  dtol = d;
              else if (m > 0)
                  dtol = tol;
              else
                  dtol = -tol;
              b = b + dtol;
              wylicz fb = f(b);
              if (!DIFFSIGN(fb, fc))
              {
                  c = a; fc = fa;
                  d = e = b - a;
              }
          }
          else break;
      }
      return b;
  }

Powyższy algorytm nie sprawdza, czy wewnątrz przedziału $\langle a,b\rangle$ występuje miejsce zerowe, zatem należy upewnić się co do tego przed jego uruchomieniem.

Natomiast, jeśli miejsc zerowych jest kilka, podawane jest tylko jedno z nich.

Przykład

Dla funkcji $f(x)={\frac {1}{x-3}}-6$ na [3,01; 4]:

• a₁ = 3,010000000000, b₁ = 4,000000000000, c₁ = 3,010000000000
• a₂ = 4,000000000000, b₂ = 3,950000000000, c₂ = 3,010000000000
• a₃ = 3,950000000000, b₃ = 3,480000000000, c₃ = 3,010000000000
• a₄ = 3,480000000000, b₄ = 3,245000000000, c₄ = 3,010000000000
• a₅ = 3,245000000000, b₅ = 3,127500000000, c₅ = 3,245000000000
• a₆ = 3,127500000000, b₆ = 3,185075000000, c₆ = 3,127500000000
• a₇ = 3,185075000000, b₇ = 3,170992625000, c₇ = 3,127500000000
• a₈ = 3,170992625000, b₈ = 3,166554383174, c₈ = 3,170992625000
• a₉ = 3,166554383174, b₉ = 3,166669581069, c₉ = 3,166554383174
• a₁₀ = 3,166669581069, b₁₀ = 3,166666668630, c₁₀ = 3,166554383174
• a₁₁ = 3,166666668630, b₁₁ = 3,166666666667, c₁₁ = 3,166666668630
• a₁₂ = 3,166666666667, b₁₂ = 3,166666666668, c₁₂ = 3,166666666667

Własności

Metoda Brenta zachowuje się znacznie lepiej niż metoda Dekkera (algorytm A), jednak przewyższa ją metoda Dekkera (algorytm R). Przykładowo dla powyższej funkcji potrzeba 12 kroków (licząc z początkowym), podczas gdy metoda Dekkera R potrzebuje tylko 5 kroków.

Co więcej, gdy przetestujemy funkcję $f(x)=(x+3)(x-1)^{2}$ mającą zera w punktach –3 i 1 z różnymi punktami początkowymi i końcowymi na przedziale $[-5;5]$ z siatką co 0.01, to dla niektórych zakresów będzie wymagane aż 69 kroków jak $[-4.8;1.8],$ podczas gdy dla metody Dekkera R żaden zakres nie będzie wymagał więcej niż 14 kroków.