Metody Numeryczne - Równanie sześcienne

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:

Jeden pierwiastek rzeczywisty, dwa pierwiastki zespolone:
Pierwiastek rzeczywisty potrójny:
Trzy pierwiastki rzeczywiste, w tym dwa podwójne:
Trzy różne pierwiastki rzeczywiste:

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

x₁, x₂, x₃

–

rozwiązania

Zmienne pomocnicze

delta

–

wyróżnik

Lista kroków

K01:	x₁ ← 0	Jeden pierwiastek zawsze jest równy 0
K02:		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

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 - x₀ dzieli bez reszty wielomian W ( x ), czyli:

Oznacza to, iż x₀ 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 = x₀. 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

Jeśli znamy pierwiastek rzeczywisty x₀ równania sześciennego ax³+bx²+cx+d = 0, to przy pomocy schematu Hornera dzielimy wielomian przez dwumian x - x₀ 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:

Przy pomocy metody Newtona znajdujemy pierwiastek rzeczywisty x₀równania sześciennego ax³+bx²+cx+d = 0.
Dzielimy wielomian ax³+bx²+cx+d przez dwumian x - x₀ za pomocą schematu Hornera i otrzymujemy wielomian a'x²+b'x+c'.
Znajdujemy pozostałe dwa pierwiastki x₁ i x₂ równania kwadratowego a'x²+b'x+c' = 0. Pierwiastki te są jednocześnie pierwiastkami równania sześciennego ax³+bx²+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

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 ) = ax³+ bx² + 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: (x_k1, f ( x_k1 ) ) oraz K2: (x_k2, f ( x_k2 ) ). 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:

Tylko jeden punkt krytyczny
Dwa różne punkty krytyczne:
Brak punktów krytycznych:

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

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:

Jeśli wartość f ( x_k ) jest różna od zera, to pierwiastek x₀ leży na lewo lub na prawo punktu x_k w zależności od znaków współczynnika a oraz f ( x_k ):

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

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

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

Obliczamy wartości funkcji f ( x_k1 ) i f ( x_k2 ). Jeśli f ( x_k1 ) = 0, to x_k1 jest pierwiastkiem rzeczywistym. Jeśli f ( x_k2 ) = 0, to x_k2 jest pierwiastkiem rzeczywistym.
Jeśli żadna z wartości funkcji f ( x_k1 ) i f ( x_k2 ) 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:

W przeciwnym razie pierwiastek x₀ leży:
Na lewo od x_k1, jeśli znaki a i f ( x_k1 ) są zgodne:

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

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

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 x₀. Jeśli nie, to wyznaczamy punkt startowy x₀ 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 x₀. Jeśli tak, to przechodzimy do kroku 3 z pierwiastkiem x₀. Jeśli nie, to sprawdzamy, czy wartości funkcji w punktach krytycznych mają różne znaki. Jeśli tak, to za punkt startowy x₀ wybieramy punkt pomiędzy punktami krytycznymi i przechodzimy do kroku 2. Inaczej punkt startowy x₀ wybieramy na lewo od x_k1 lub na prawo od x_k2 w zależności od znaków współczynnika a oraz wartości funkcji w punkcie x_k1 i przechodzimy do kroku 2.
Jeśli nie ma punktów krytycznych, to za punkt startowy x₀ wybieramy 0 i przechodzimy do kroku 2.
Metodą Newtona szukamy przybliżonego pierwiastka równania sześciennego, przyjmując za punkt startowy x₀ ustawione w kroku 1. Wynik umieszczamy w x₀_.
Mając pierwiastek x₀ dzielimy wielomian sześcienny równania przez dwumian x - x₀ za pomocą schematu Hornera.
Rozwiązujemy wynikowe równanie kwadratowe, otrzymując pozostałe dwa pierwiastki x₁ i x₂.

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

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:

Trzy pierwiastki rzeczywiste, 2 podwójne:
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 - x₁ i rozwiązując otrzymane w wyniku dzielenia równanie kwadratowe, jak robiliśmy to już w algorytmach 2 i 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

Równania

Równanie sześcienne

Algorytm nr 1

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

Dane wejściowe:

Dane wyjściowe

Zmienne pomocnicze

Lista kroków

Dzielenie wielomianów schematem Hornera

Algorytm nr 2

Algorytm nr 3

Algorytm nr 4