Metoda Dekkera

Pomysłem Theodorusa Dekkera z 1969 roku było połączenie metody równego podziału (bisekcji) z metodą siecznych. Przypuśćmy, że chcemy rozwiązać równanie f(x) = 0. Podobnie jak w metodzie bisekcji musimy zainicjować metodę Dekkera dwoma punktami, powiedzmy x0 i x1, takimi że f(x0) i f(x1) mają różne znaki. Jeśli funkcja f jest ciągła na [x0,x1] wtedy twierdzenie Darboux gwarantuje że zero leży pomiędzy x0 i x1.

Opis algorytmu A[edytuj | edytuj kod]

W każdej iteracji zaangażowane są trzy punkty:

b – ostatnie przybliżenie wyniku,

c – kontrapunkt b, czyli taki punkt w którym funkcja f ma przeciwny znak niż w punkcie b,

a – poprzednia wartość punktu a, używany jest on do wyliczania następnego punktu metodą siecznych.

Metoda bisekcji natomiast używa punktów b i c tworząc punkt m pomiędzy nimi na środku przedziału. Wyliczany jest ciąg xi, którego ostatni element oznaczany jest przez x, a poprzedni przez xp. Oprócz tego jest punkt xk, który jest ostatnim punktem w ciągu, który ma różny znak niż x. Następny punkt x wyliczany jest dwoma metodami: siecznych i bisekcji i wybierany jest ten obliczony z metody siecznych jeśli leży pomiędzy punktem b (ze względów dokładnościowych z pewną poprawką) a punktem m wyliczonym z bisekcji.

Następnie badanie czy f(x) czy f(xk) leży bliżej zera i jeśli f(x) leży bliżej zera wtedy b ma wartość x a c xk, w przeciwnym razie zamiana.

Algorytm A[edytuj | edytuj kod]

\delta (b)

= |b*MACHINE_EPSILON| gdzie MACHINE_EPSILON jest to dokładność dla liczb zmiennoprzecinkowych, dla double jest to 2.2204460492503131e-16^[1]

DIFFSIGN(x,y)

x<=0

AND

y>=0

OR

x>=0

AND

y<=0

– sprawdzenie czy x i y są różnych znaków; sprawdzenie znaku iloczynu może dawać złe wyniki gdy iloczyn równy zeru a liczby nie; dla poniższych algorytmów lepiej tu stosować nieostre znaki < i >, tak że w przypadku gdy jedno lub dwa są zerami zostaną uznane za liczby różnych znaków.

 bool between(x, a, b)
 {
        if (b > a)
               return x >= a && x <= b;
        else
               return x >= b && x <= a;
 }
 lfun(b, a, fb, fa)
 {
        if (fb != fa)
                return b - fb*(b - a) / (fb - fa); //[[metoda siecznych]]
        else if (fa != 0)
                return INFINITY;
        else return b;
 }
 hfun(b, c)
 {
        return b + SIGN(c-b)*<math>\delta(b)</math>;
 }
 mfun(b, c)
 {
        return 0.5*(b + c); //[[Metoda równego podziału|metoda bisekcji]]
 }
 vfun(l, b, c)
 {
        h = hfun(b, c);
        m = mfun(b, c);
        if (between(l, h, m))
                return l;
        else if (|l - b| <= <math>\delta(b)</math>)
                return h;
        else
                return m;
 }
 DekkerA(x0, x1, Eps)
 {
        wylicz fxp = f(x0);
        wylicz fx = f(x1);
        if (|fx| <= |fxp|)
        {
                b = x1;
                a = c = x0;
                fa = fxp;
                fb = fx;
 }
        else
        {
                b = x0;
                a = c = x1;
                fa = fx;
                fb = fxp;
 }
        xk = x0;
        fxk = fxp;
        x = x1;
        while (|b - c| > 2 * Eps)
        {
                lambda = lfun(b, a, fb, fa);
                xp = x;
                x = vfun(lambda, b, c);
                fxp = fx;
                wylicz fx = f(x);
                if (DIFFSIGN(fxp,fx))
                {
                        xk = xp;
                        fxk = fxp;
                }
                if (|fx| <= |fxk|)
                {
                        a = b;
                        b = x; c = xk;
                        fa = fb;
                        fb = fx;
                }
                else
                {
                        b = xk; a = c = x;
                        fa = fx;
                        fb = fxk;
                }
        }
        return b;
 }

Przykład[edytuj | edytuj kod]

Weźmy funkcję $f(x)={\frac {1}{x-3}}-6$ ^[2]. Osiąga ona ∞ w punkcie 3.0 i zero w $x=3{\frac {1}{6}}.$ Testowy program w Visual C++ radzi sobie z zakresem [3;4] gdzie ∞ traktuje jako liczbę dodatnią, ale na wszelki wypadek zawęźmy przedział do [3.01; 4]. Stosujemy dokładność double i przerywamy po osiągnięciu dokładności (różnicy między b i c) 1e-12.

na początku a₁ = 3.01, b₁ = 4, f(a₁) = 94, f(b₁) = -5; obliczamy c₁ = 3.01
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.166188864569, c₈ = 3.170992625000
a₉ = 3.166188864569, b₉ = 3.166679068378, c₉ = 3.166188864569
a₁₀ = 3.166679068378, b₁₀ = 3.166666702220, c₁₀ = 3.166188864569
a₁₁ = 3.166666702220, b₁₁ = 3.166666666664, c₁₁ = 3.166666702220
a₁₂ = 3.166666666664, b₁₂ = 3.166666666667, c₁₂ = 3.166666702220
a₁₃ = 3.166666666667, b₁₃ = 3.166666666667, c₁₃ = 3.166666666667

Własności[edytuj | edytuj kod]

O ile dla powyższej funkcji mającej jeden pierwiastek z przedziale metoda zachowuje się w miarę dobrze, to na przykład dla funkcji $f(x)=(x+3)(x-1)^{2}$ mającej zera w punktach -3 i 1, szukając na przedziale [-4, 4/3] występują kłopoty: b zbliża się powoli jednostronnie do zera, a c ma cały czas wartość -4 jeszcze w 74-tej pętli, aż dopiero w następnej (dla double) b osiąga wartość dokładnie 1, f(b) osiąga zero, i c staje się równe b. Innym problemem jest przypadek, gdy w pobliżu pierwiastka pochodne dążą do zera, co powoduje bardzo powolną zbieżność. Brent podaje taką funkcję: $f(x)=x\cdot \ exp(-x^{-2})$ jeśli x różne od zera i zero dla zera.

Opis algorytmu M[edytuj | edytuj kod]

W tej metodzie bada się to kiedy ostatnio przedział |b-c| został zmniejszony o połowę. Tutaj ustawiana jest zmienna age. Dodatkowo, gdy age=3, wtedy wołana jest funkcja rfun, która umożliwia zbieżność szybciej niż lfun.

Algorytm M[edytuj | edytuj kod]

between, lfun, hfun, mfun, $\delta (b)$ oraz DIFFSIGN zdefiniowane jak w algorytmie A.

  wfun(l, b, c)
  {
      h = hfun(b, c);
      m = mfun(b, c);
      if (between(l, h, m))
          return l;
      else if (f|l - b| <= <math>\delta(b)</math> && !between(l, b, m))
          return h;
      else
          return m;
  }
  ffun(a, b, fa, fb)
  {
      return (fa - fb) / (a - b);
  }
  rfun(b, a, d, fb, fa, fd)
  {
      alpha = ffun(b, d, fb, fd)*fa;
      beta = ffun(a, d, fa, fd)*fb;
      if (beta != alpha)
          return b - beta*(b - a) / (beta - alpha);
      else if (alpha != 0)
          return INFINITY;
      else
          return 0; //beta == alpha == 0
  }
  DekkerM(x0, x1, Eps)
  {
      itercnt = 0;
      wylicz fxp = f(x0);
      wylicz fx = f(x1);
      if (x0 == x1) return fx;
      if (!DIFFSIGN(fx, fxp)) return 0;
      if (|fx| <= |fxp|)
      {
          b = x1;
          a = c = x0;
          fa = fxp;
          fb = fx;
      }
      else
      {
          b = x0;
          a = c = x1;
          fa = fx;
          fb = fxp;
      }
      xk = xp = x0;
      fxk = fxp;
      x = x1;
      itercnt = 1;
      int age = 0;
      bp = b;
      cp = c;
      ap = a;
      while (|b - c| > 2 * Eps)
      {
          itercnt++;
          age++;
          if (|b - c| <= (0.5 + 2 * MACHINE_EPSILON)*(|bp - cp| + <math>\delta(b)</math>) age = 1;
          xp = x;
          if (age<=2)
          {
              lambda = lfun(b, a, fb, fa);
              if (|lambda - b| < <math>\delta(b)</math>) break;
              x = wfun(lambda, b, c);
          }
          else if (age == 3)
          {
              rho = rfun(b, a, d, fb, fa, fd);
              if (|rho - b| < <math>\delta(b)</math>) break;
              x = wfun(rho, b, c);
          }
          else
              x = mfun(b, c);
          fxp = fx;
          wylicz fx = f(x);
          if (DIFFSIGN(fxp,fx))
          {
              xk = xp;
              fxk = fxp;
          }
          bp = b;
          fbp = fb;
          ap = a;
          fap = fa;
          cp = c;
          if (|fx| <= |fxk|)
          {
              a = b;
              fa = fb;
              b = x;
              fb = fx;
              c = xk;
          }
          else
          {
              b = xk;
              fb = fxk;
              a = c = x;
              fa = fx;
          }
          if (b == x || b == bp)
          {
              d = ap;
              fd = fap;
          }
          else
          {
              d = bp;
              fd = fbp;
          }
      }
      return b;
  }

Przykład[edytuj | edytuj kod]

Weźmy funkcję $f(x)=(x+3)(x-1)^{2}$ mającej zera w punktach -3 i 1, będziemy szukać na przedziale [-4, 4/3]

a₀ = -4.000000000000, b₀ = 1.333333333333, c₀ = -4.000000000000
age = 1 lfun a₂ = 1.333333333333, b₂ = 1.232558139535, c₂ = -4.000000000000
age = 2 lfun a₃ = 1.232558139535, b₃ = 1.141223295850, c₃ = -4.000000000000
age = 3 rfun a₄ = 1.141223295850, b₄ = 1.070756096437, c₄ = -4.000000000000
age = 4 bisekcja a₅ = 1.070756096437, b₅ = -1.464621951782, c₅ = -4.000000000000
age = 1 lfun a₆ = -1.464621951782, b₆ = -2.732310975891, c₆ = -4.000000000000
age = 1 lfun a₇ = -3.366155487945, b₇ = -2.732310975891, c₇ = -3.366155487945
age = 1 lfun a₈ = -2.732310975891, b₈ = -2.953018236685, c₈ = -3.366155487945
age = 2 lfun a₉ = -2.953018236685, b₉ = -3.007123150382, c₉ = -2.953018236685
age = 1 lfun a₁₀ = -3.007123150382, b₁₀ = -2.999830139829, c₁₀ = -3.007123150382
age = 1 lfun a₁₁ = -2.999830139829, b₁₁ = -2.999999396604, c₁₁ = -3.007123150382
age = 2 lfun a₁₂ = -2.999999396604, b₁₂ = -3.000000000051, c₁₂ = -2.999999396604
age = 1 lfun a₁₃ = -3.000000000051, b₁₃ = -3.000000000000, c₁₃ = -3.000000000051

Dla $f(x)={\frac {1}{x-3}}-6$ na [3,01; 4] potrzeba 12 kroków.

Opis algorytmu R[edytuj | edytuj kod]

Częściej wołana jest funkcja rfun zamiast lfun co dla wielu funkcji powoduje znaczny przyrost zbieżności.

Algorytm R[edytuj | edytuj kod]

Identyczny jak Algorytm M z wyjątkiem fragmentu wyliczania x:

    ..........................................................
    if (itercnt == 2)
    {
        lambda = lfun(b, a, fb, fa);
        if (|lambda - b| < <math>\delta(b)</math>) break;
        x = wfun(lambda, b, c);
    }
    else if (itercnt >= 3 && age <= 3)
    {
        rho = rfun(b, a, d, fb, fa, fd);
        if (|rho - b| < <math>\delta(b)</math>) break;
        x = wfun(rho, b, c);
    }
    else if (itercnt >= 3 && age == 4)
    {
        rho = rfun(b, a, d, fb, fa, fd);
        if (|rho - b| < <math>\delta(b)</math>) break;
        x = wfun(2 * rho - b, b, c);
    }
    else
        x = mfun(b, c);
    ..........................................................

Wikibooks:pl:Kody źródłowe/Metoda Dekkera