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:

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:

C++

// Równanie sześcienne
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne 0042
//---------------------------

#include <iostream>
#include <windows.h>
#include <cmath>
#include <iomanip>

using namespace std;

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

// Wyświetla + przed
// liczbą dodatnią
// '-' przed ujemną
//------------------
void pl(double x)
{
  if(x >= 0)
    cout << '+';
  else
    cout << '-';
  cout << fabs(x);
}

// 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 << "\nRównanie: "
       << a << "*x^3";
  pl(b);
  cout << "*x^2";
  pl(c);
  cout << "*x = 0\n";

  // Jeden z pierwiastków
  // zawsze jest równy 0
  cout << "x1 = +0.0000\n";

  // Rozwiązujemy równanie
  // kwadratowe
  delta = b * b - 4 * a * c;

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

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

int main()
{
  SetConsoleOutputCP(CP_UTF8);
  SetConsoleCP(CP_UTF8);

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

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

  cout << endl;
  system("pause");
  return 0;
}

Wynik

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

Równanie: 3.000*x^3+12.000*x^2+15.000*x = 0
x1 = +0.000
x2 = -2.000-1.000i
x3 = -2.000+1.000i

Równanie: 9.000*x^3+12.000*x^2+4.000*x = 0
x1 = +0.000
x2 = -0.667
x3 = -0.667

Równanie: 1.000*x^3-1.000*x^2-2.000*x = 0
x1 = +0.000
x2 = -1.000
x3 = +2.000

Python (dodatek)

# Równanie sześcienne
# (C)2019 mgr Jerzy Wałaszek
# Metody numeryczne 0042
#---------------------------

import math

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

# Funkcja wypisuje liczbę
# ze znakiem + lub -
#------------------------
def pl(x):
    if x >= 0:
        print('+', end='')
    else:
        print('-', end='')
    print(f"{math.fabs(x):.3f}",
          end='')

# Funkcja rozwiązująca
# równanie sześcienne
# a,b,c - współczynniki
#         wielomianu
#----------------------
def ce(a, b, c):
    print()
    print(f"Równanie: {a:.4f}*x^3",
                end='')
    pl(b)
    print("*x^2", end='')
    pl(c)
    print("*x = 0")

    # Jeden z pierwiastków
    # zawsze jest równy 0
    print("x1 = +0.0000")

    # Rozwiązujemy równanie
    # kwadratowe
    delta = b * b - 4 * a * c

    if abs(delta) <= eps:
        c_val = - b / 2 / a
        print("x2 = ", end='')
        pl(c_val)
        print("\nx3 = ", end='')
        pl(c_val)
    elif delta < -eps:
        re = - b / 2 / a
        im = math.sqrt(-delta) / 2 / a
        print("x2 = ", end='')
        pl(re)
        pl(-im)
        print("i\nx3 = ", end='')
        pl(re)
        pl(im)
        print("i", end='')
    else:
        delta = math.sqrt(delta)
        print("x2 = ", end='')
        pl((-b - delta) / 2 / a)
        print("\nx3 = ", end='')
        pl((-b + delta) / 2 / a)
    print()

# Program główny
#---------------

print("Rozwiązywanie równań "
      "sześciennych typu\n"
      "  3    2\n"
      "ax + bx + cx = 0\n"
      "----------------")

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

print()
input("Naciśnij Enter...")

do podrozdziału do 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:

do podrozdziału do 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.

C++

// Równanie sześcienne
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne 0043
//---------------------------

#include <iostream>
#include <windows.h>
#include <cmath>
#include <iomanip>

using namespace std;

// Tutaj definiujemy
// dane wejściowe
//------------------

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

// Wyświetla + przed
// liczbą dodatnią
// '-' przed ujemną
//------------------
void pl(double x)
{
  if(x >= 0)
    cout << '+';
  else
    cout << '-';
  cout << fabs(x);
}

// 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 << "\nRównanie: ";
  pl(a);
  cout << "*x^3";
  pl(b);
  cout << "*x^2";
  pl(c);
  cout << "*x";
  pl(d);
  cout << " = 0\n";

  // 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)
    {
      // Zaznaczamy, że wynik
      // jest zespolony
      result = false;
      // Część rzeczywista
      x1 = -bp/2/ap;
      // Część urojona
      x2 = sqrt(-delta)/2/ap;
    }
    else
    {
      delta = sqrt(delta);
      x1 = (-bp-delta)/2/ap;
      x2 = (-bp+delta)/2/ap;
    }

    // Wyświetlamy wyniki
    cout << "x0 = ";
    pl(x0);
    cout << endl;
    if(result)
    {
      cout << "x1 = ";
      pl(x1);
      cout << "\nx2 = ";
      pl(x2);
    }
    else
    {
      cout << "x1 = ";
      pl(x1);
      pl(-x2);
      cout << "i\n"
              "x2 = ";
      pl(x1);
      pl(x2);
      cout << "i";
    }
  }
  else
    cout << "Błąd, przekroczono "
            "zadaną liczbę obiegów\n"
            "OBLICZENIA PRZERWANE!";
  cout << endl;
}

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

int main()
{
  SetConsoleOutputCP(CP_UTF8);
  SetConsoleCP(CP_UTF8);

  cout << setprecision(3)
       << fixed;
  cout << "Równania sześcienne typu\n"
          "  3    2\n"
          "ax + bx + cx + d = 0\n"
          "--------------------\n";

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

  cout << endl;
  system("pause");
  return 0;
}

Wynik

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

Równanie: +1.000*x^3-6.000*x^2+11.000*x-6.000 = 0
x0 = +1.000
x1 = +2.000
x2 = +3.000

Równanie: +1.000*x^3-1.000*x^2-16.000*x-20.000 = 0
x0 = -2.000
x1 = -2.000
x2 = +5.000

Równanie: +1.000*x^3-2.000*x^2-10.000*x-25.000 = 0
x0 = +5.000
x1 = -1.500-1.658i
x2 = -1.500+1.658i

Python (dodatek)

# Równanie sześcienne
# (C)2026 mgr Jerzy Wałaszek
# Metody numeryczne 0043
# ---------------------------

import math

# Tutaj definiujemy
# dane wejściowe
# ------------------

# Dokładność wyznaczania
# pierwiastka
epsx = 1e-14
# Dokładność porównania
# z zerem
epsy = 1e-14

# Funkcja wypisuje liczbę
# ze znakiem + lub -
#------------------------
def pl(x):
    if x >= 0:
        print('+', end='')
    else:
        print('-', end='')
    print(f"{math.fabs(x):.3f}",
          end='')

# Funkcja rozwiązująca
# równanie sześcienne
# a,b,c,d - współczynniki
#           wielomianu
# ------------------------
def nh(a, b, c, d):
    n = 32
    result = False

    print("\nRównanie: ", end="")
    pl(a)
    print("*x^3", end="")
    pl(b)
    print("*x^2", end="")
    pl(c)
    print("*x", end="")
    pl(d)
    print(" = 0")

    # Najpierw szukamy pierwiastka
    # rzeczywistego metodą Newtona
    x0 = 1  # Punkt startowy
    while n > 1:
        n -= 1
        # Obliczamy wartość funkcji
        # w punkcie x0
        f0 = x0*(c+x0*(b+x0*a))+ d

        # Sprawdzamy, czy funkcja
        # jest dostatecznie
        # bliska zeru
        if abs(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 abs(x1 - x0) < epsx:
            result = True
            break

        # Kontynuujemy obliczenia

    if result:
        # Mamy pierwiastek x0.
        # Schematem Hornera
        # dzielimy wielomian
        # równania przez
        # dwumian (x - x0)

        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 abs(delta) < epsx:
            x1 = -bp / 2 / ap
            x2 = x1
        elif delta < 0:
            # Zaznaczamy, że wynik
            # jest zespolony
            result = False
            # Część rzeczywista
            x1 = -bp / 2 / ap
            # Część urojona
            x2 = math.sqrt(-delta)/2/ap
        else:
            delta = math.sqrt(delta)
            x1 = (-bp-delta)/2/ap
            x2 = (-bp+delta)/2/ap

        # Wyświetlamy wyniki
        print("x0 = ", end="")
        pl(x0)
        print()
        if result:
            print("x1 = ", end="")
            pl(x1)
            print("\nx2 = ", end="")
            pl(x2)
        else:
            print("x1 = ", end="")
            pl(x1)
            pl(-x2)
            print("i\nx2 = ", end="")
            pl(x1)
            pl(x2)
            print("i")
    else:
        print("Błąd, przekroczono "
              "zadaną liczbę obiegów\n"
              "OBLICZENIA PRZERWANE!")
    print()

# Program główny
# ---------------

print("Równania sześcienne typu\n"
      "  3    2\n"
      "ax + bx + cx + d = 0\n"
      "--------------------")

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

input("Naciśnij Enter...")

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.000*x^3+0.000*x^2+0.000*x+0.000 = 0
x0 = +0.000
x1 = -0.000-0.000i
x2 = -0.000+0.000i


Równanie: +1.000*x^3+3.000*x^2+3.000*x+1.000 = 0
x0 = -1.000
x1 = -1.000-0.000i
x2 = -1.000+0.000i


Równanie: +1.000*x^3-6.000*x^2+12.000*x-8.000 = 0
x0 = +2.000
x1 = +2.000-0.000i
x2 = +2.000+0.000i


Równanie: +1.000*x^3+9.000*x^2+27.000*x+27.000 = 0
Błąd, przekroczono zadaną liczbę obiegów
OBLICZENIA PRZERWANE!

do podrozdziału do 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

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 algorytmu:

C++

// Równanie sześcienne
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne 0044
//---------------------------

#include <iostream>
#include <windows.h>
#include <cmath>
#include <iomanip>

using namespace std;

// Tutaj definiujemy
// dane wejściowe
//------------------

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

// Wyświetla + przed
// liczbą dodatnią
// '-' przed ujemną
//------------------
void pl(double x)
{
  if(x >= 0)
    cout << '+';
  else
    cout << '-';
  cout << fabs(x);
}

// 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 << "\nRównanie: ";
  pl(a);
  cout << "*x^3";
  pl(b);
  cout << "*x^2";
  pl(c);
  cout << "*x";
  pl(d);
  cout << " = 0\n";

  // 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)
      // Sygnalizujemy
      // pierwiastek x0
      result = true;
    else if(a*f0 > 0)
      // Punkt startowy
      // na lewo od x0
      x0 -= 1;
    else
      // Punkt startowy
      // na prawo od x0
      x0 += 1;
  }
  else if(delta < 0)
    // Brak punktów krytycznych,
    // dowolny punkt startowy
    x0 = 0;
  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)
      // Sygnalizujemy
      // pierwiastek x0
      result = true;
    else if(fabs(f1) < epsy)
    {
      x0 = x1;
      // Sygnalizujemy
      // pierwiastek x0
      result = true;
    }
    else if(f0*f1 < 0)
      x0 = (x0+x1)/2;
    else
    {
      if(a*f0 > 0)
        // Punkt startowy
        // na lewo od x0
        x0 -= 1;
      else
        // Punkt startowy
        // na prawo od x1
        x0 = x1+1;
    }
  }
  // Jeśli nie ma pierwiastka,
  // to szukamy go metodą Newtona
  if(!result)
    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)
    {
      // Sygnalizujemy
      // wynik zespolony
      result = false;
      // Część rzeczywista
      x1 = -bp/2/ap;
      // Część urojona
      x2 = sqrt(-delta)/2/ap;
    }
    else
    {
      delta = sqrt(delta);
      x1 = (-bp-delta)/2/ap;
      x2 = (-bp+delta)/2/ap;
    }

    // Wyświetlamy wyniki
    cout << "x0 = ";
    pl(x0);
    cout << endl;
    if(result)
    {
      cout << "x1 = ";
      pl(x1);
      cout << "\nx2 = ";
      pl(x2);
    }
    else
    {
      cout << "x1 = ";
      pl(x1);
      pl(-x2);
      cout << "i\n"
              "x2 = ";
      pl(x1);
      pl(x2);
      cout << "i";
    }
  }
  else
    cout << "Błąd, przekroczono "
            "zadaną liczbę obiegów\n"
            "OBLICZENIA PRZERWANE!";
  cout << endl;
}

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

int main()
{
  SetConsoleOutputCP(CP_UTF8);
  SetConsoleCP(CP_UTF8);

  cout << setprecision(3)
       << fixed;
  cout <<
  "Równania sześcienne typu\n"
  "  3    2\n"
  "ax + bx + cx + d = 0\n"
  "--------------------\n";

  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;
  system("pause");
  return 0;
}

Wynik

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

Równanie: +1.000*x^3+0.000*x^2+0.000*x+0.000 = 0
x0 = +0.000
x1 = +0.000
x2 = +0.000

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

Równanie: +1.000*x^3-6.000*x^2+12.000*x-8.000 = 0
x0 = +2.000
x1 = +2.000
x2 = +2.000

Równanie: +1.000*x^3+9.000*x^2+27.000*x+27.000 = 0
x0 = -3.000
x1 = -3.000
x2 = -3.000

Równanie: +1.000*x^3-6.000*x^2+11.000*x-6.000 = 0
x0 = +2.000
x1 = +1.000
x2 = +3.000

Równanie: +1.000*x^3-1.000*x^2-16.000*x-20.000 = 0
x0 = -2.000
x1 = -2.000
x2 = +5.000

Równanie: +1.000*x^3-2.000*x^2-10.000*x-25.000 = 0
x0 = +5.000
x1 = -1.500-1.658i
x2 = -1.500+1.658i

Python (dodatek)

# Równanie sześcienne
# (C)2019 mgr Jerzy Wałaszek
# Metody numeryczne 0044
#---------------------------

import math

# Tutaj definiujemy
# dane wejściowe
#------------------

# Dokładność wyznaczania
# pierwiastka
epsx = 1e-14
# Dokładność porównania
# z zerem
epsy = 1e-14

# Funkcja wypisuje liczbę
# ze znakiem + lub -
#------------------------
def pl(x):
    if x >= 0:
        print('+', end='')
    else:
        print('-', end='')
    print(f"{math.fabs(x):.3f}",
          end='')

# Funkcja rozwiązująca
# równanie sześcienne
# a,b,c,d - współczynniki
#       wielomianu
#------------------------
def nh(a, b, c, d):
    n = 32
    result = False

    print("\nRównanie: ", end="")
    pl(a)
    print("*x^3", end="")
    pl(b)
    print("*x^2", end="")
    pl(c)
    print("*x", end="")
    pl(d)
    print(" = 0")

    # Znajdujemy punkty
    # krytyczne
    delta = b*b-3*a*c
    if abs(delta) < epsx:
        x0 = -b/3/a
        f0 = x0*(c+x0*(b+x0*a))+d
        if abs(f0) < epsy:
            # Sygnalizujemy
            # pierwiastek x0
            result = True
        elif a * f0 > 0:
            # Punkt startowy
            # na lewo od x0
            x0 -= 1
        else:
            # Punkt startowy
            # na prawo od x0
            x0 += 1
    elif delta < 0:
        # Brak punktów
        # krytycznych,
        # dowolny punkt
        # startowy
        x0 = 0.0
    else:
        delta = math.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 abs(f0) < epsy:
            # Sygnalizujemy
            # pierwiastek x0
            result = True
        elif abs(f1) < epsy:
            x0 = x1
            # Sygnalizujemy
            # pierwiastek x0
            result = True
        elif f0 * f1 < 0:
            x0 = (x0+x1)/2
        else:
            if a * f0 > 0:
                # Punkt startowy
                # na lewo od x0
                x0 -= 1
            else:
                # Punkt startowy
                # na prawo od x1
                x0 = x1+1

    # Jeśli nie ma
    # pierwiastka, to
    # szukamy go metodą
    # Newtona
    if not result:
        while n > 0:
            n -= 1
            # Obliczamy wartość
            # funkcji w punkcie x0
            f0 = x0*(c+x0*(b+x0*a))+d

            # Sprawdzamy, czy funkcja
            # jest dostatecznie
            # bliska zeru
            if abs(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 abs(x1 - x0) < epsx:
                result = True
                break

            # Kontynuujemy
            # obliczenia

    if result:
        # Mamy pierwiastek
        # x0. Schematem
        # Hornera dzielimy
        # wielomian równania
        # przez dwumian
        # (x - x0)
        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 abs(delta) < epsx:
            x1 = -bp/2/ ap
            x2 = x1
        elif delta < 0:
            # Sygnalizujemy
            # wynik zespolony
            result = False
            # Część rzeczywista
            x1 = -bp/2/ap
            # Część urojona
            x2 = math.sqrt(-delta)/2/ap
        else:
            delta = math.sqrt(delta)
            x1 = (-bp-delta)/2/ap
            x2 = (-bp+delta)/2/ap

        # Wyświetlamy wyniki
        print("x0 = ", end="")
        pl(x0)
        print()
        if result:
            print("x1 = ", end="")
            pl(x1)
            print("\nx2 = ", end="")
            pl(x2)
            print()
        else:
            print("x1 = ", end="")
            pl(x1)
            pl(-x2)
            print("i\nx2 = ", end="")
            pl(x1)
            pl(x2)
            print("i")
    else:
        print("Błąd, przekroczono\n"
              "zadaną liczbę obiegów\n"
              "OBLICZENIA PRZERWANE!")

# Program główny
#---------------

print("Równania sześcienne typu\n"
      "  3    2\n"
      "ax + bx + cx + d = 0\n"
      "--------------------")

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)

print()
input("Naciśnij Enter...")

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.

do podrozdziału do 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:

C++

// Równanie sześcienne
// (C)2019 mgr Jerzy Wałaszek
// Metody numeryczne 045
//---------------------------

#include <iostream>
#include <windows.h>
#define _USE_MATH_DEFINES
#include <cmath>
#include <iomanip>

using namespace std;

// Tutaj definiujemy
// dane wejściowe
//------------------

// Dokładność porównania
// z zerem
double epsy = 1e-12;

// Wyświetla + przed
// liczbą dodatnią,
// '-' przed ujemną
//------------------
void pl(double x)
{
  if(x >= 0)
    cout << '+';
  else
    cout << '-';
  cout << fabs(x);
}

// 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;
  // Znacznik wyniku
  // zespolonego
  bool cmplx = false;

  cout << "\nRównanie: ";
  pl(a);
  cout << "*x^3";
  pl(b);
  cout << "*x^2";
  pl(c);
  cout << "*x";
  pl(d);
  cout << " = 0\n";

  // 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.0;
  q = 2*a*a*a/27.0-a*b/3.0+c;
  delta = p*p*p/27.0+q*q/4.0;

  // Sprawdzamy przypadki

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

    // Schematem Hornera
    // dzielimy wielomian
    // równania sześciennego
    // przez dwumian (x - x1)
    double aa,bb,cc;
    aa = 1.0;
    bb = a+aa*x1;
    cc = b+bb* x1;

    // Obliczamy wyróżnik
    // równania kwadratowego
    delta = bb*bb-4.0*aa*cc;
    delta = sqrt(-delta);
    // Część rzeczywista
    x2 = -bb/2.0/aa;
    // Część urojona
    x3 = delta/2.0/aa;
  }
  else
  { // 3 pierwiastki rzeczywiste
    a /= 3.0;
    p = sqrt(-p);
    double pp = 2.0/sqrt(3.0)*p;
    q = 3.0*sqrt(3.0)*q/2.0/p/p/p;
    x1 =  pp*sin(asin(q)/3.0)-a;
    x2 = -pp*sin(asin(q)/3.0+M_PI/3.0)-a;
    x3 =  pp*cos(asin(q)/3.0+M_PI/6.0)-a;
  }
  cout << "x1 = ";
  pl(x1);
  cout << endl;
  if(cmplx)
  {
    cout << "x2 = ";
    pl(x2);
    pl(-x3);
    cout << "i"
            "\nx3 = ";
    pl(x2);
    pl(x3);
    cout << "i";
  }
  else
  {
    cout << "x2 = ";
    pl(x2);
    cout << "\nx3 = ";
    pl(x3);
  }
  cout << endl;
}

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

int main()
{
  SetConsoleOutputCP(CP_UTF8);
  SetConsoleCP(CP_UTF8);

  cout << setprecision(3)
       << fixed;
  cout << "Równania sześcienne typu\n"
          "  3    2\n"
          "ax + bx + cx + d = 0\n"
          "--------------------\n";

  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;
  system("pause");
  return 0;
}

Wynik

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

Równanie: +1.000*x^3+0.000*x^2+0.000*x+0.000 = 0
x1 = +0.000
x2 = +0.000
x3 = +0.000

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

Równanie: +1.000*x^3-6.000*x^2+12.000*x-8.000 = 0
x1 = +2.000
x2 = +2.000
x3 = +2.000

Równanie: +1.000*x^3+9.000*x^2+27.000*x+27.000 = 0
x1 = -3.000
x2 = -3.000
x3 = -3.000

Równanie: +1.000*x^3-6.000*x^2+11.000*x-6.000 = 0
x1 = +2.000
x2 = +1.000
x3 = +3.000

Równanie: +1.000*x^3-1.000*x^2-16.000*x-20.000 = 0
x1 = +5.000
x2 = -2.000
x3 = -2.000

Równanie: +1.000*x^3-2.000*x^2-10.000*x-25.000 = 0
x1 = +5.000
x2 = -1.500-1.658i
x3 = -1.500+1.658i

Python (dodatek)

# Równanie sześcienne
# (C)2026 mgr Jerzy Wałaszek
# Metody numeryczne 045
#---------------------------

import math

# Tutaj definiujemy
# dane wejściowe
#------------------

# Dokładność porównania
# z zerem
epsy = 1e-12

# Stała PI
PI = math.pi

# Wyświetla + przed
# liczbą dodatnią,
# '-' przed ujemną
#------------------
def pl(x):
    if x >= 0:
        print("+", end="")
    else:
        print("-", end="")
    print(f"{abs(x):.3f}", end="")

# Funkcja rozwiązująca
# równanie sześcienne
# a,b,c,d - współczynniki
#       wielomianu
#------------------------
def ce(a, b, c, d):
    cmplx = False

    print("\nRównanie: ", end="")
    pl(a)
    print("*x^3", end="")
    pl(b)
    print("*x^2", end="")
    pl(c)
    print("*x", end="")
    pl(d)
    print(" = 0")

    # 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 abs(delta) < epsy:
        # 3 pierwiastki,
        # 2 podwójne
        sgn = (q > 0)-(q < 0)
        q = sgn*(abs(q)/2)**(1/3)
        x1 = -2*q-a/3
        x2 = q-a/3
        x3 = x2
    elif delta > 0:
        # 1 pierwiastek
        # rzeczywisty,
        # 2 zespolone
        cmplx = True
        q /= -2
        delta = math.sqrt(delta)
        v1 = q + delta
        s1 = (v1 > 0)-(v1 < 0)
        v2 = q - delta
        s2 = (v2 > 0)-(v2 < 0)
        x1 = (s1*abs(v1)**(1/3)+
              s2*abs(v2)**(1/3)-
              a/3)

        # Schematem Hornera
        # dzielimy wielomian
        # równania
        # sześciennego
        # przez dwumian
        # (x - x1)
        aa = 1.0
        bb = a+aa*x1
        cc = b+bb*x1

        # Obliczamy
        # wyróżnik
        # równania
        # kwadratowego
        delta = bb*bb-4*aa*cc
        delta = math.sqrt(-delta)
        # Część rzeczywista
        x2 = -bb/2/aa
        # Część urojona
        x3 = delta/2/aa
    else:
        # 3 pierwiastki
        # rzeczywiste
        a /= 3.0
        p = math.sqrt(-p)
        pp = 2/math.sqrt(3)*p
        q = (3*math.sqrt(3)
            *q/2/p/p/p)
        q = max(-1,min(1,q))
        ang = math.asin(q)/3
        x1 =  pp*math.sin(ang)-a
        x2 = -pp*math.sin(ang+PI/3)-a
        x3 =  pp*math.cos(ang+PI/6)-a

    print("x1 = ", end="")
    pl(x1)
    print()
    if cmplx:
        print("x2 = ", end="")
        pl(x2)
        pl(-x3)
        print("i\nx3 = ", end="")
        pl(x2)
        pl(x3)
        print("i")
    else:
        print("x2 = ", end="")
        pl(x2)
        print("\nx3 = ", end="")
        pl(x3)
        print()

# Program główny
#---------------

print("Równania sześcienne typu\n"
      "  3    2\n"
      "ax + bx + cx + d = 0\n"
      "--------------------")

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)

print()
input("Naciśnij Enter...")

do podrozdziału do strony

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.

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