Serwis Edukacyjny
w I-LO w Tarnowie
obrazek

Materiały dla uczniów liceum

  Wyjście       Spis treści       Wstecz       Dalej  

obrazek

Autor artykułu: mgr Jerzy Wałaszek

©2019 mgr Jerzy Wałaszek
I LO w Tarnowie

obrazek

Równania

Równanie sześcienne

SPIS TREŚCI
Podrozdziały

Algorytm nr 1

Funkcja sześcienna (ang. qubic function) jest wielomianem trzeciego stopnia:

Równanie sześcienne (ang. qubic equation) powstaje, gdy przyrównamy funkcję sześcienną do 0:

Współczynniki a, b, c i d są liczbami rzeczywistymi.

Wykres funkcji sześciennej wygląda następująco:

obrazek

W zależności od położenia tego wykresu mamy następujące możliwości:

  1. Jeden pierwiastek rzeczywisty, dwa pierwiastki zespolone:
        obrazek
     
  2. Pierwiastek rzeczywisty potrójny:
        obrazek
     
  3. Trzy pierwiastki rzeczywiste, w tym dwa podwójne:
        obrazek
     
  4. Trzy różne pierwiastki rzeczywiste:
        obrazek

Równanie sześcienne zawsze posiada przynajmniej jeden pierwiastek rzeczywisty. Zrozumiesz to natychmiast po przyjrzeniu się wykresowi funkcji sześciennej.

Rozwiązanie równań sześciennych jest dużo trudniejsze obliczeniowo od rozwiązywania równań kwadratowych. Dokładne metody opracowano dopiero w XVI wieku. Zajmijmy się najpierw prostymi przypadkami.

Mamy równanie sześcienne typu:

W równaniu tym współczynnik d jest równy zero. Skoro tak, to zmienną x możemy wyprowadzić przed nawias:

Wynika z tego, iż jeden z pierwiastków tego równania jest równy 0. W nawiasie otrzymaliśmy wielomian kwadratowy, który rozwiązujemy algorytmem z poprzedniego rozdziału. W ten sposób otrzymamy dwa dalsze pierwiastki, które mogą być rzeczywiste lub zespolone.

Algorytm rozwiązywania równania sześciennego przy d=0

Dane wejściowe:

ε dokładność przyrównania do zera
a,b,c współczynniki, a różne od 0

Dane wyjściowe

x1, x2, x3 rozwiązania

Zmienne pomocnicze

delta wyróżnik

Lista kroków

K01: x1 ← 0 Jeden pierwiastek zawsze jest równy 0
K02: obrazek Obliczamy wyróżnik delta
K03: Jeśli |delta| < ε, to

Zakończ
delta = 0? Pierwiastek podwójny
K04: Jeśli delta < 0, to

Zakończ
delta < 0? Dwa pierwiastki zespolone
K05: delta > 0. Dwa pierwiastki rzeczywiste
K06: Zakończ  

Poniższy program rozwiązuje równania sześcienne przy założeniu, iż współczynnik d jest równy 0, a współczynnik a jest różny od zera:

Przykładowy program w języku C++
// Równanie sześcienne
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne
//---------------------------

#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

// Tutaj definiujemy dane wejściowe
//----------------------------------
double eps = 1e-12;  // Dokładność porównania z zerem

// Funkcja rozwiązująca równanie sześcienne
// a,b,c - współczynniki wielomianu
//-----------------------------------------
void ce(double a, double b, double c)
{
    double delta;

    cout << endl;
    cout << "Równanie: " << a << " * x^3 + " << b << " * x^2 + " << c << " * x = 0" << endl;

    // Jeden z pierwiastków zawsze jest równy 0

    cout << "x1 = 0" << endl;

    // Rozwiązujemy równanie kwadratowe

    delta = b * b - 4 * a * c;

    if(fabs(delta) < eps) cout << "x2 = x3 = " << - b / 2 / a;
    else if(delta < 0)
    {
        double re = - b / 2 / a;
        double im = sqrt(-delta) / 2 / a;
        cout << "x2 = " << re << " + " << -im << "i" << endl
             << "x3 = " << re << " + " << im << "i";
    }
    else
    {
        delta = sqrt(delta);
        cout << "x2 = " << (-b - delta) / 2 / a << endl
             << "x3 = " << (-b + delta) / 2 / a;
    }
    cout << endl;
}

// Program główny
//---------------

int main()
{
    setlocale(LC_ALL,"");

    cout << fixed;
    cout << "Rozwiązywanie równań sześciennych typu" << endl
         << "  3    2" << endl
         << "ax + bx + cx = 0" << endl
         << "----------------" << endl;

    ce(3,12,15);
    ce(9,12,4);
    ce(1,-1,-2);

    cout << endl;

    return 0;
}
Wynik
Rozwiązywanie równań sześciennych typu
  3    2
ax + bx + cx = 0
----------------

Równanie: 3.000000 * x^3 + 12.000000 * x^2 + 15.000000 * x = 0
x1 = 0
x2 = -2.000000 + -1.000000i
x3 = -2.000000 + 1.000000i

Równanie: 9.000000 * x^3 + 12.000000 * x^2 + 4.000000 * x = 0
x1 = 0
x2 = x3 = -0.666667

Równanie: 1.000000 * x^3 + -1.000000 * x^2 + -2.000000 * x = 0
x1 = 0
x2 = -1.000000
x3 = 2.000000
Na początek:  podrozdziału   strony 

Dzielenie wielomianów schematem Hornera

Wielomiany można dzielić, podobnie jak liczby. Wynikiem dzielenia wielomianu W ( x ) przez wielomian X ( x ) jest wielomian P ( x ) oraz wielomian reszty z dzielenia R ( x ) o następujących własnościach:

Jeśli wielomian reszty z dzielenia R( x ) = 0, to wielomian W( x ) jest podzielny przez wielomian X ( x ), a wynikiem dzielenia jest P ( x ).

Załóżmy, że

W takim przypadku:

Problem sprowadza się do wyznaczenia współczynników a', b', c' oraz d'. Współczynniki te w prosty sposób wyznaczamy wg tzw. schematu Hornera:

Jeśli d' = 0, to dwumian x - x0 dzieli bez reszty wielomian W ( x ), czyli:

Oznacza to, iż x0 jest pierwiastkiem wielomianu W ( x ).

Jeśli d' jest różne od 0, to możemy potraktować je jako wartość wielomianu W ( x ) dla x = x0. W tej interpretacji schemat Hornera może służyć do szybkiego wyznaczania wartości wielomianu (korzystamy z tego w programach).

Przykład:

Na początek:  podrozdziału   strony 

Algorytm nr 2

Jeśli znamy pierwiastek rzeczywisty x0 równania sześciennego ax3+bx2+cx+d = 0, to przy pomocy schematu Hornera dzielimy wielomian przez dwumian x - x0 i pozostałe dwa pierwiastki wyznaczamy algorytmem rozwiązywania równań kwadratowych. Ponieważ każde równanie sześcienne ma co najmniej jeden pierwiastek rzeczywisty, to do jego znalezienia możemy wykorzystać jedną z metod znajdowania przybliżonego pierwiastka funkcji. Szczególnie dobrym wyborem będzie tutaj metoda Newtona, ponieważ jest szybka i pochodną wielomianu liczy się bardzo prosto:

Algorytm będzie zatem następujący:

  1. Przy pomocy metody Newtona znajdujemy pierwiastek rzeczywisty x0 równania sześciennego ax3+bx2+cx+d = 0.
  2. Dzielimy wielomian ax3+bx2+cx+d przez dwumian x - x0 za pomocą schematu Hornera i otrzymujemy wielomian a'x2+b'x+c'.
  3. Znajdujemy pozostałe dwa pierwiastki x1 i x2 równania kwadratowego a'x2+b'x+c' = 0. Pierwiastki te są jednocześnie pierwiastkami równania sześciennego ax3+bx2+cx+d = 0.

Poniższy program jest przykładem zastosowania tego algorytmu.

Przykładowy program w języku C++
// Równanie sześcienne
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne
//---------------------------

#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

// Tutaj definiujemy dane wejściowe
//----------------------------------
double epsx = 1e-14;  // Dokładność wyznaczania pierwiastka
double epsy = 1e-14;  // Dokładność porównania z zerem

// Funkcja rozwiązująca równanie sześcienne
// a,b,c,d - współczynniki wielomianu
//-----------------------------------------
void nh(double a, double b, double c, double d)
{
    double x0,x1,x2,f0,f1;

    int n = 32;

    bool result = false;

    cout << endl
         << "Równanie: "
         << a << " * x^3 + "
         << b << " * x^2 + "
         << c << " * x + "
         << d << " = 0" << endl;

    // Najpierw szukamy pierwiastka rzeczywistego metodą Newtona

    x0 = 1; // Punkt startowy

    while(--n)
    {
       // Obliczamy wartość funkcji w punkcie x0

       f0 = x0 * (c + x0 * (b + x0 * a)) + d;

       // Sprawdzamy, czy funkcja jest dostatecznie bliska zeru

       if(fabs(f0) < epsy)
       {
           result = true;
           break;
       }

       // Obliczamy wartość pierwszej pochodnej funkcji

       f1 = x0 * (2 * b + 3 * a * x0) + c;

       // Zapamiętujemy bieżące przybliżenie

       x1 = x0;

       // Obliczamy kolejne przybliżenie

       x0 -= f0/f1;

       // Sprawdzamy, czy odległość pomiędzy dwoma ostatnimi przybliżeniami
       // jest mniejsza od założonej dokładności

       if(fabs(x1 - x0) < epsx)
       {
           result = true;
           break;
       }

       // Kontynuujemy obliczenia
    }

    if(result)
    {
        // Mamy pierwiastek x0
        // Schematem Hornera dzielimy wielomian równania
        // przez dwumian (x - x0)

        double ap,bp,cp,delta;

        ap = a;
        bp = b + ap * x0;
        cp = c + bp * x0;

        // Rozwiązujemy równanie ap * x^2 + bp * x + cp = 0

        delta = bp * bp - 4 * ap * cp;

        if(fabs(delta) < epsx)
        {
            x1 = x2 = - bp / 2 / ap;
        }
        else if(delta < 0)
        {
            result = false; // Zaznaczamy, że wynik jest zespolony
            x1 = -bp / 2 / ap; // Część rzeczywista
            x2 = sqrt(-delta) / 2 / ap; // Część urojona
        }
        else
        {
            delta = sqrt(delta);
            x1 = (-bp - delta) / 2 / ap;
            x2 = (-bp + delta) / 2 / ap;
        }

        // Wyświetlamy wyniki

        cout << "x0 = " << x0 << endl;
        if(result) cout << "x1 = " << x1 << endl
                        << "x2 = " << x2;
        else       cout << "x1 = " << x1 << " + " << -x2 << "i" << endl
                        << "x2 = " << x1 << " + " << x2 << "i";
    }
    else cout << "Błąd, przekroczono zadaną liczbę obiegów" << endl
              << "OBLICZENIA PRZERWANE!";
    cout << endl;
}

// Program główny
//---------------

int main()
{
    setlocale(LC_ALL,"");

    cout << fixed;
    cout << "Równania sześcienne typu" << endl
         << "  3    2" << endl
         << "ax + bx + cx + d = 0" << endl
         << "--------------------" << endl;

    nh(1,-6,11,-6);
    nh(1,-1,-16,-20);
    nh(1,-2,-10,-25);

    cout << endl;

    return 0;
}
Wynik
Równania sześcienne typu
  3    2
ax + bx + cx + d = 0
--------------------

Równanie: 1.000000 * x^3 + -6.000000 * x^2 + 11.000000 * x + -6.000000 = 0
x0 = 1.000000
x1 = 2.000000
x2 = 3.000000

Równanie: 1.000000 * x^3 + -1.000000 * x^2 + -16.000000 * x + -20.000000 = 0
x0 = -2.000000
x1 = -2.000000
x2 = 5.000000

Równanie: 1.000000 * x^3 + -2.000000 * x^2 + -10.000000 * x + -25.000000 = 0
x0 = 5.000000
x1 = -1.500000 + -1.658312i
x2 = -1.500000 + 1.658312i

Zastosowana metoda Newtona może zawieść dla niektórych równań sześciennych, ponieważ w pewnych wypadkach pochodna zeruje się (metodę Newtona stosuje się wtedy, gdy pierwsza pochodna jest różna od zera pomiędzy punktem startowym, a najbliższym pierwiastkiem), co wprowadza duże błędy do wyniku. Na przykład, dla równań:

otrzymujemy wyniki:

Równania sześcienne typu
  3    2
ax + bx + cx + d = 0
--------------------

Równanie: 1.000000 * x^3 + 0.000000 * x^2 + 0.000000 * x + 0.000000 = 0
x0 = 0.000018
x1 = -0.000009 + -0.000015i
x2 = -0.000009 + 0.000015i

Równanie: 1.000000 * x^3 + 3.000000 * x^2 + 3.000000 * x + 1.000000 = 0
x0 = -0.999984
x1 = -1.000008 + -0.000014i
x2 = -1.000008 + 0.000014i

Równanie: 1.000000 * x^3 + -6.000000 * x^2 + 12.000000 * x + -8.000000 = 0
x0 = 1.999982
x1 = 2.000009 + -0.000015i
x2 = 2.000009 + 0.000015i

Równanie: 1.000000 * x^3 + 9.000000 * x^2 + 27.000000 * x + 27.000000 = 0
x0 = -2.999979
x1 = -3.000010 + -0.000018i
x2 = -3.000010 + 0.000018i
Na początek:  podrozdziału   strony 

Algorytm nr 3

Kolejny algorytm jest ulepszeniem algorytmu nr 2, w którym pozbywamy się z rachunków pochodnej równej zero.  Przyjrzyjmy się wykresowi funkcji sześciennej  f ( x ) = ax3 + bx2 + cx + d :

obrazek

Zwróć uwagę, iż funkcja f ( x ) posiada dwa punkty, w których styczne do wykresu (zaznaczone kolorem czerwonym) są równoległe do osi OX, a zatem nie przecinają jej lub leżą na tej osi. Punkty te oznaczyliśmy symbolami K1 i K2 i będziemy nazywać punktami krytycznymi (ang. critical points) funkcji f ( x ). Mają one współrzędne K1: (xk1, f ( xk1 )) oraz K2: (xk2, f ( xk2 )). Styczna do wykresu funkcji w punkcie K1 lub K2 będzie równoległa do osi OX, jeśli pierwsza pochodna funkcji f(x) ma w tym punkcie wartość zero. Ta własność pozwala nam znaleźć punkty K1 i K2:

Otrzymujemy równanie kwadratowe. Wyliczamy wyróżnik delta:

W zależności od wartości wyróżnika delta mamy trzy przypadki:

  1. Tylko jeden punkt krytyczny

     
  2. Dwa różne punkty krytyczne:

     
  3. Brak punktów krytycznych:

Teraz rozważymy kolejno możliwe położenia pierwiastków równania sześciennego.

  1. Jeden punkt krytyczny:

    Sprawdzamy, czy funkcja ma w tym punkcie wartość 0. Jeśli tak, to punkt ten jest potrójnym pierwiastkiem równania sześciennego:
    obrazek
    Jeśli wartość  f ( xk ) jest różna od zera, to pierwiastek x0 leży na lewo lub na prawo punktu xk w zależności od znaków współczynnika a oraz  f ( xk ):

    Na lewo od xk: znaki współczynnika a i  f ( xk ) są zgodne:
    obrazek

    Na prawo od xk: znaki współczynnika a i f(xk) są przeciwne:
    obrazek

    Informację o położeniu pierwiastka x0 wykorzystujemy do określenia punktu startowego dla metody Newtona.
     
  2. Dwa punkty krytyczne:

    Obliczamy wartości funkcji  f ( xk1 ) i  f ( xk2 ). Jeśli  f ( xk1 ) = 0, to xk1 jest pierwiastkiem rzeczywistym. Jeśli  f ( xk2 ) = 0, to xk2 jest pierwiastkiem rzeczywistym.
    Jeśli żadna z wartości funkcji  f ( xk1 )f ( xk2 ) nie jest równa 0, to sprawdzamy, czy wartości te są różnych znaków. Jeśli tak, to pierwiastek leży pomiędzy nimi:
    obrazek
    W przeciwnym razie pierwiastek x0 leży:
    Na lewo od xk1, jeśli znaki af ( xk1 ) są zgodne:
    obrazek

    Na prawo od xk2, jeśli znaki a i  f ( xk1 ) są przeciwne:
    obrazek
  3. Brak punktów krytycznych:
    Punkt startowy wybieramy dowolnie.

Z powyższych rozważań wyłania się następujący schemat algorytmu:

  1. Wyznaczamy punkty krytyczne:
    Jeśli jest tylko jeden punkt krytyczny, to sprawdzamy, czy jest on pierwiastkiem. Jeśli tak, to przechodzimy do kroku 3 z pierwiastkiem x0. Jeśli nie, to wyznaczamy punkt startowy x0 na lewo lub na prawo od punktu krytycznego wg znaków współczynnika a oraz wartości funkcji w punkcie krytycznym i przechodzimy do kroku 2.
    Jeśli są dwa punkty krytyczne, to sprawdzamy, czy jeden z nich jest pierwiastkiem x0. Jeśli tak, to przechodzimy do kroku 3 z pierwiastkiem x0. Jeśli nie, to sprawdzamy, czy wartości funkcji w punktach krytycznych mają różne znaki. Jeśli tak, to za punkt startowy x0 wybieramy punkt pomiędzy punktami krytycznymi i przechodzimy do kroku 2. Inaczej punkt startowy x0 wybieramy na lewo od xk1 lub na prawo od xk2 w zależności od znaków współczynnika a oraz wartości funkcji w punkcie xk1 i przechodzimy do kroku 2.
    Jeśli nie ma punktów krytycznych, to za punkt startowy x0 wybieramy 0 i przechodzimy do kroku 2.
  2. Metodą Newtona szukamy przybliżonego pierwiastka równania sześciennego, przyjmując za punkt startowy x0 ustawione w kroku 1. Wynik umieszczamy w x0.
  3. Mając pierwiastek x0 dzielimy wielomian sześcienny równania przez dwumian x - x0 za pomocą schematu Hornera.
  4. Rozwiązujemy wynikowe równanie kwadratowe, otrzymując pozostałe dwa pierwiastki x1 i x2.

Poniższy program jest przykładową realizacją przedstawionego powyżej algorytmu:

Przykładowy program w języku C++
// Równanie sześcienne
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne
//---------------------------

#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

// Tutaj definiujemy dane wejściowe
//----------------------------------
double epsx = 1e-14;  // Dokładność wyznaczania pierwiastka
double epsy = 1e-14;  // Dokładność porównania z zerem

// Funkcja rozwiązująca równanie sześcienne
// a,b,c,d - współczynniki wielomianu
//-----------------------------------------
void nh(double a, double b, double c, double d)
{
    double x0,x1,x2,f0,f1,delta;

    int n = 32;

    bool result = false;

    cout << endl
         << "Równanie: "
         << a << " * x^3 + "
         << b << " * x^2 + "
         << c << " * x + "
         << d << " = 0" << endl;

    // Znajdujemy punkty krytyczne

    delta = b * b - 3 * a * c;

    if(fabs(delta) < epsx)
    {
        x0 = - b / 3 / a;
        f0 = x0 * (c + x0 * (b + x0 * a)) + d;
        if(fabs(f0) < epsy) result = true; // Sygnalizujemy pierwiastek x0
        else if(a * f0 > 0) x0 -= 1;  // Punkt startowy na lewo od x0
        else                x0 += 1;  // Punkt startowy na prawo od x0
    }
    else if(delta < 0) x0 = 0; // Brak punktów krytycznych, dowolny punkt startowy
    else
    {
        delta = sqrt(delta);
        x0 = (-b - delta) / 3 / a;
        x1 = (-b + delta) / 3 / a;
        f0 = x0 * (c + x0 * (b + x0 * a)) + d;
        f1 = x1 * (c + x1 * (b + x1 * a)) + d;
        if(fabs(f0) < epsy) result = true; // Sygnalizujemy pierwiastek x0
        else if(fabs(f1) < epsy)
        {
           x0 = x1;
           result = true;                  // Sygnalizujemy pierwiastek x0
        }
        else if(f0 * f1 < 0) x0 = (x0 + x1) / 2;
        else
        {
            if(a * f0 > 0) x0 -= 1;     // Punkt startowy na lewo od x0
            else           x0 = x1 + 1; // Punkt startowy na prawo od x1
        }
    }
    if(!result)  // Jeśli nie ma pierwiastka, to szukamy go metodą Newtona
        while(--n)
        {
           // Obliczamy wartość funkcji w punkcie x0

           f0 = x0 * (c + x0 * (b + x0 * a)) + d;

           // Sprawdzamy, czy funkcja jest dostatecznie bliska zeru

           if(fabs(f0) < epsy)
           {
               result = true;
               break;
           }

           // Obliczamy wartość pierwszej pochodnej funkcji

           f1 = x0 * (2 * b + 3 * a * x0) + c;

           // Zapamiętujemy bieżące przybliżenie

           x1 = x0;

           // Obliczamy kolejne przybliżenie

           x0 -= f0/f1;

           // Sprawdzamy, czy odległość pomiędzy dwoma ostatnimi przybliżeniami
           // jest mniejsza od założonej dokładności

           if(fabs(x1 - x0) < epsx)
           {
               result = true;
               break;
           }

           // Kontynuujemy obliczenia
        }

    if(result)
    {
        // Mamy pierwiastek x0
        // Schematem Hornera dzielimy wielomian równania
        // przez dwumian (x - x0)

        double ap,bp,cp;

        ap = a;
        bp = b + ap * x0;
        cp = c + bp * x0;

        // Rozwiązujemy równanie ap * x^2 + bp * x + cp = 0

        delta = bp * bp - 4 * ap * cp;

        if(fabs(delta) < epsx) x1 = x2 = - bp / 2 / ap;
        else if(delta < 0)
        {
            result = false;             // Sygnalizujemy wynik zespolony
            x1 = -bp / 2 / ap;          // Część rzeczywista
            x2 = sqrt(-delta) / 2 / ap; // Część urojona
        }
        else
        {
            delta = sqrt(delta);
            x1 = (-bp - delta) / 2 / ap;
            x2 = (-bp + delta) / 2 / ap;
        }

        // Wyświetlamy wyniki

        cout << "x0 = " << x0 << endl;
        if(result) cout << "x1 = " << x1 << endl
                        << "x2 = " << x2;
        else       cout << "x1 = " << x1 << " + " << -x2 << "i" << endl
                        << "x2 = " << x1 << " + " << x2 << "i";
    }
    else cout << "Błąd, przekroczono zadaną liczbę obiegów" << endl
              << "OBLICZENIA PRZERWANE!";
    cout << endl;
}

// Program główny
//---------------

int main()
{
    setlocale(LC_ALL,"");

    cout << fixed;
    cout << "Równania sześcienne typu" << endl
         << "  3    2" << endl
         << "ax + bx + cx + d = 0" << endl
         << "--------------------" << endl;

    nh(1,0,0,0);
    nh(1,3,3,1);
    nh(1,-6,12,-8);
    nh(1,9,27,27);
    nh(1,-6,11,-6);
    nh(1,-1,-16,-20);
    nh(1,-2,-10,-25);

    cout << endl;

    return 0;
}
Wynik
Równania sześcienne typu
  3    2
ax + bx + cx + d = 0
--------------------

Równanie: 1.000000 * x^3 + 0.000000 * x^2 + 0.000000 * x + 0.000000 = 0
x0 = -0.000000
x1 = -0.000000
x2 = -0.000000

Równanie: 1.000000 * x^3 + 3.000000 * x^2 + 3.000000 * x + 1.000000 = 0
x0 = -1.000000
x1 = -1.000000
x2 = -1.000000

Równanie: 1.000000 * x^3 + -6.000000 * x^2 + 12.000000 * x + -8.000000 = 0
x0 = 2.000000
x1 = 2.000000
x2 = 2.000000

Równanie: 1.000000 * x^3 + 9.000000 * x^2 + 27.000000 * x + 27.000000 = 0
x0 = -3.000000
x1 = -3.000000
x2 = -3.000000

Równanie: 1.000000 * x^3 + -6.000000 * x^2 + 11.000000 * x + -6.000000 = 0
x0 = 2.000000
x1 = 1.000000
x2 = 3.000000

Równanie: 1.000000 * x^3 + -1.000000 * x^2 + -16.000000 * x + -20.000000 = 0
x0 = -2.000000
x1 = -2.000000
x2 = 5.000000

Równanie: 1.000000 * x^3 + -2.000000 * x^2 + -10.000000 * x + -25.000000 = 0
x0 = 5.000000
x1 = -1.500000 + -1.658312i
x2 = -1.500000 + 1.658312i

Porównaj otrzymane wyniki z wynikami poprzedniego programu. Zwróć uwagę, iż teraz pierwiastki potrójne są znajdowane prawidłowo. Pierwiastki niekoniecznie otrzymujemy w kolejności od najmniejszego do największego. Jeśli ma to znaczenie, to należy zastosować prosty algorytm sortujący.

Na początek:  podrozdziału   strony 

Algorytm nr 4

Równania sześcienne można rozwiązywać, podobnie jak równania kwadratowe, licząc wyróżnik delta. Jednakże w tym przypadku postać wyróżnika jest dużo bardziej skomplikowana. Nie będziemy wyprowadzać wzorów, opiszemy jedynie metodę rozwiązywania dowolnych równań sześciennych.

Mamy równanie sześcienne:

Dzielimy wszystkie współczynniki przez a:

Przyjmujemy nowe współczynniki:

i otrzymujemy równanie wynikowe:

W pierwszym kroku wyliczamy wyróżnik delta. Ponieważ ma on dosyć skomplikowany wzór, zastosujemy podstawienia:

Wyróżnik równania zdefiniujemy teraz następująco:

W zależności od wartości wyróżnika delta mamy trzy możliwe przypadki:


  1. Trzy pierwiastki rzeczywiste, 2 podwójne:

     

  2. Tylko jeden pierwiastek rzeczywisty:

    Pozostałe dwa pierwiastki zespolone znajdujemy, jeśli to konieczne, dzieląc schematem Hornera wielomian równania sześciennego przez dwumian x - x1 i rozwiązując otrzymane w wyniku dzielenia równanie kwadratowe, jak robiliśmy to już w algorytmach 2 i 3.
     

  3. Trzy pierwiastki rzeczywiste:

Poniższy program jest przykładową realizacją przedstawionego powyżej algorytmu:

Przykładowy program w języku C++
// Równanie sześcienne
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne
//---------------------------

#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

// Tutaj definiujemy dane wejściowe
//----------------------------------
double epsy = 1e-12;  // Dokładność porównania z zerem
const double M_PI = 3.14159265358979323846;

// Funkcja rozwiązująca równanie sześcienne
// a,b,c,d - współczynniki wielomianu
//-----------------------------------------
void ce(double a, double b, double c, double d)
{
    double x1,x2,x3,p,q,delta;
    bool cmplx = false;     // Znacznik wyniku zespolonego

    cout << endl
         << "Równanie: "
         << a << " * x^3 + "
         << b << " * x^2 + "
         << c << " * x + "
         << d << " = 0" << endl;

    // Modyfikujemy współczynniki równania

    d /= a;
    c /= a;
    b /= a;
    a = b;
    b = c;
    c = d;

    // Obliczamy wyróżnik delta

    p = b - a * a / 3;
    q = 2 * a * a * a / 27 - a * b / 3 + c;
    delta = p * p * p / 27 + q * q / 4;

    // Sprawdzamy przypadki

    if(fabs(delta) < epsy) // 3 pierwiastki, 2 podwójne
    {
        q = ((q > 0) - (q < 0)) * pow(fabs(q) / 2, 1 / 3.0);
        x1 = -2 * q - a / 3;
        x2 = x3 = q - a / 3;
    }
    else if(delta > 0)    // 1 pierwiastek rzeczywisty, 2 zespolone
    {
        cmplx = true;     // Wynik zespolony
        q /= -2;
        delta = sqrt(delta);
        x1 = ((q + delta > 0) - (q + delta < 0)) * pow(fabs(q + delta), 1 / 3.0) +
             ((q - delta > 0) - (q - delta < 0)) * pow(fabs(q - delta), 1 / 3.0) - a / 3;

        // Schematem Hornera dzielimy wielomian równania sześciennego
        // przez dwumian x - x1

        double aa,bb,cc;

        aa = 1;
        bb = a + aa * x1;
        cc = b + bb * x1;

        // Obliczamy wyróżnik równania kwadratowego

        delta = bb * bb - 4 * aa * cc;

        delta = sqrt(-delta);

        x2 = - bb / 2 / aa ; // Część rzeczywista
        x3 = delta / 2 / aa; // Część urojona
    }
    else  // 3 pierwiastki rzeczywiste
    {
        a /= 3;
        p = sqrt(-p);
        double pp = 2 / sqrt(3) * p;
        q = 3 * sqrt(3) * q / 2 / p / p / p;
        x1 = pp * sin(asin(q) / 3) - a;
        x2 = - pp * sin(asin(q) / 3 + M_PI / 3) - a;
        x3 = pp * cos(asin(q) / 3 + M_PI / 6) - a;
    }
    cout << "x1 = " << x1 << endl;
    if(cmplx) cout << "x2 = " << x2 << " + " << -x3 << "i" << endl
                   << "x3 = " << x2 << " + " << x3 << "i";
    else      cout << "x2 = " << x2 << endl
                   << "x3 = " << x3;
    cout << endl;
}

// Program główny
//---------------

int main()
{
    setlocale(LC_ALL,"");

    cout << fixed;
    cout << "Równania sześcienne typu" << endl
         << "  3    2" << endl
         << "ax + bx + cx + d = 0" << endl
         << "--------------------" << endl;

    ce(1,0,0,0);
    ce(1,3,3,1);
    ce(1,-6,12,-8);
    ce(1,9,27,27);
    ce(1,-6,11,-6);
    ce(1,-1,-16,-20);
    ce(1,-2,-10,-25);

    cout << endl;

    return 0;
}
Wynik
Równania sześcienne typu
  3    2
ax + bx + cx + d = 0
--------------------

Równanie: 1.000000 * x^3 + 0.000000 * x^2 + 0.000000 * x + 0.000000 = 0
x1 = -0.000000
x2 = 0.000000
x3 = 0.000000

Równanie: 1.000000 * x^3 + 3.000000 * x^2 + 3.000000 * x + 1.000000 = 0
x1 = -1.000000
x2 = -1.000000
x3 = -1.000000

Równanie: 1.000000 * x^3 + -6.000000 * x^2 + 12.000000 * x + -8.000000 = 0
x1 = 2.000000
x2 = 2.000000
x3 = 2.000000

Równanie: 1.000000 * x^3 + 9.000000 * x^2 + 27.000000 * x + 27.000000 = 0
x1 = -3.000000
x2 = -3.000000
x3 = -3.000000

Równanie: 1.000000 * x^3 + -6.000000 * x^2 + 11.000000 * x + -6.000000 = 0
x1 = 2.000000
x2 = 1.000000
x3 = 3.000000

Równanie: 1.000000 * x^3 + -1.000000 * x^2 + -16.000000 * x + -20.000000 = 0
x1 = 5.000000
x2 = -2.000000
x3 = -2.000000

Równanie: 1.000000 * x^3 + -2.000000 * x^2 + -10.000000 * x + -25.000000 = 0
x1 = 5.000000
x2 = -1.500000 + -1.658312i
x3 = -1.500000 + 1.658312i
Na początek:  podrozdziału   strony 

Zespół Przedmiotowy
Chemii-Fizyki-Informatyki

w I Liceum Ogólnokształcącym
im. Kazimierza Brodzińskiego
w Tarnowie
ul. Piłsudskiego 4
©2019 mgr Jerzy Wałaszek

Materiały tylko do użytku dydaktycznego. Ich kopiowanie i powielanie jest dozwolone
pod warunkiem podania źródła oraz niepobierania za to pieniędzy.

Pytania proszę przesyłać na adres email: i-lo@eduinf.waw.pl

Serwis wykorzystuje pliki cookies. Jeśli nie chcesz ich otrzymywać, zablokuj je w swojej przeglądarce.
Informacje dodatkowe.