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 czwartego stopnia

SPIS TREŚCI
Podrozdziały

Algorytm nr 1

Równanie czwartego stopnia (ang. quartic equation) zbudowane jest z wielomianu czwartego stopnia:

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

Im wyższy stopień wielomianu, tym bardziej skomplikowane wzory opisują jego rozwiązania. Równanie stopnia czwartego jest ostatnim, które można rozwiązać algebraicznie. Udowodniono, że równania powyżej stopnia czwartego nie posiadają ogólnej metody rozwiązywania, co nie oznacza braku pierwiastków. Równania takie, oprócz przypadków szczególnych, rozwiązuje się przy pomocy metod przybliżonych.

Rozwiązania równań czwartego stopnia matematycy opracowali dopiero w XVI wieku. Powstało wiele metod rozwiązywania tych równań. Tutaj opiszemy tylko kilka najbardziej znanych.

Na początek zajmiemy się szczególnymi przypadkami równania czwartego stopnia, które stosunkowo łatwo rozwiązać (wykorzystując poznane wcześniej metody rozwiązywania równań wielomianowych).

Równanie dwukwadratowe (ang. biquadratic equation) ma postać:

Dokonujemy podstawienia:

Otrzymaliśmy równanie kwadratowe ze zmienną y. Rozwiązujemy je i otrzymujemy dwa pierwiastki: y1 i y2. Następnie otrzymujemy pierwiastki równania dwukwadratowego:

Jeśli któryś z pierwiastków y1 lub y2 jest ujemny lub zespolony, to otrzymamy w wyniku zespolone pierwiastki xi. Jeśli interesują nas tylko rozwiązania rzeczywiste, pomijamy rozwiązania zespolone. Możemy również wykonać rachunki na liczbach zespolonych i pomijać części urojone równe zero, otrzymując w ten sposób pierwiastki rzeczywiste (każda liczba rzeczywista ma część urojoną równą zero).

Załóżmy, że mamy liczbę zespoloną o postaci:

Pierwiastek zespolony tej liczby obliczymy ze wzoru:

Część rzeczywista pierwiastka zespolonego jest zawsze nieujemna. Znak części urojonej jest taki sam jak znak części urojonej pierwiastkowanej liczby (funkcja sgn(x) ma wartość -1 dla argumentów ujemnych, 0 dla argumentu 0 i 1 dla argumentu dodatniego).

Powyższy wzór pozwala nam rozwiązać dowolne równanie dwukwadratowe za pomocą metody rozwiązywania równań kwadratowych bez rozważania różnych przypadków.

Przykładowy program rozwiązuje dowolne równanie dwukwadratowe.

Przykładowy program w języku C++
// Równanie czwartego stopnia
// (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

// Definicja typu zespolonego
struct cmplx
{
    double re, im;
};

// Funkcja wylicza pierwiastek zespolony
//--------------------------------------
void csqrt(cmplx y, cmplx & x)
{
    double a,b;
    a = y.re;
    b = y.im;
    double c = sqrt(a * a + b * b);
    x.re = sqrt((c + a) / 2);
    x.im = sqrt((c - a) / 2);
    if(b < 0) x.im = - x.im;
}

// Funkcja wyświetlająca liczbę zespoloną
//---------------------------------------
void cmplxprint(cmplx x)
{
    cout << x.re; // Część rzeczywistą wyświetlamy zawsze
    if(fabs(x.im) > epsy) // Jeśli część urojona jest różna od zera,
        cout << " + " << x.im << "i"; // to wyświetlamy ją
    cout << endl;
}

// Funkcja rozwiązująca równanie dwukwadratowe
// a,b,c - współczynniki wielomianu
//-----------------------------------------
void be(double a, double b, double c)
{
    cout << endl
         << "Równanie: "
         << a << " * x^4 + "
         << b << " * x^2 + "
         << c << " = 0" << endl;

    double delta;
    cmplx y1,y2,x1,x2,x3,x4;

    // Zerujemy części urojone zmiennych

    y1.im = y2.im = x1.im = x2.im = x3.im = x4.im = 0;

    // Obliczamy wyróżnik

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

    // Jeśli wyróżnik jest nieujemny, wyliczamy pierwiastki rzeczywiste

    if(fabs(delta) < epsy) y1.re = y2.re = -b / 2 / a;
    else if(delta > 0)
    {
        delta = sqrt(delta);
        y1.re = (-b - delta) / 2 / a;
        y2.re = (-b + delta) / 2 / a;
    }
    else // Inaczej wyliczamy pierwiastki zespolone
    {
        delta = sqrt(-delta);
        y1.re = y2.re = -b / 2 / a;
        y1.im = -delta / 2 / a;
        y2.im =  delta / 2 / a;
    }

    // Obliczamy pierwiastki zespolone rozwiązań równania kwadratowego
    csqrt(y1,x1); x1.re = -x1.re; x1.im = - x1.im; // Ten jest ujemny
    csqrt(y1,x2);
    csqrt(y2,x3); x3.re = -x3.re; x3.im = - x3.im; // Ten jest ujemny
    csqrt(y2,x4);

    // Wyświetlamy wyniki:
    cout << "x1 = "; cmplxprint(x1);
    cout << "x2 = "; cmplxprint(x2);
    cout << "x3 = "; cmplxprint(x3);
    cout << "x4 = "; cmplxprint(x4);
}

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

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

    cout << fixed;
    cout << "Równania dwukwadratowe typu" << endl
         << "  4    2" << endl
         << "ax + bx + c = 0" << endl
         << "---------------------------" << endl;

    be(1,-5,4);
    be(1,-10,9);
    be(4,7,-2);
    be(3,2,5);

    cout << endl;

    return 0;
}
Wynik
Równania dwukwadratowe typu
  4    2
ax + bx + c = 0
---------------------------

Równanie: 1.000000 * x^4 + -5.000000 * x^2 + 4.000000 = 0
x1 = -1.000000
x2 = 1.000000
x3 = -2.000000
x4 = 2.000000

Równanie: 1.000000 * x^4 + -10.000000 * x^2 + 9.000000 = 0
x1 = -1.000000
x2 = 1.000000
x3 = -3.000000
x4 = 3.000000

Równanie: 4.000000 * x^4 + 7.000000 * x^2 + -2.000000 = 0
x1 = -0.000000 + -1.414214i
x2 = 0.000000 + 1.414214i
x3 = -0.500000
x4 = 0.500000

Równanie: 3.000000 * x^4 + 2.000000 * x^2 + 5.000000 = 0
x1 = -0.691976 + 0.901201i
x2 = 0.691976 + -0.901201i
x3 = -0.691976 + -0.901201i
x4 = 0.691976 + 0.901201i
Na początek:  podrozdziału   strony 

Algorytm nr 2

Równanie czwartego stopnia o postaci:

możemy przekształcić do postaci:

Równanie to posiada pierwiastek równy zero. Pozostałe pierwiastki są pierwiastkami równania sześciennego i znajdujemy je metodami opisanymi w poprzednim rozdziale.

Przykładowy program w języku C++
// Równanie czwartego stopnia
// (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 czwartego stopnia
// a,b,c,d - współczynniki wielomianu
//-----------------------------------------
void qe(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^4 + "
         << b << " * x^3 + "
         << c << " * x^2 + "
         << d << " * x = 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 = 2 * 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 = 0" << endl << "x2 = " << x1 << endl;
    if(cmplx) cout << "x3 = " << x2 << " + " << -x3 << "i" << endl
                   << "x4 = " << x2 << " + " << x3 << "i";
    else      cout << "x3 = " << x2 << endl
                   << "x4 = " << x3;
    cout << endl;
}

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

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

    cout << fixed;
    cout << "Równania czwartego stopnia typu" << endl
         << "  4    3    2" << endl
         << "ax + bx + cx + dx = 0" << endl
         << "-------------------------------" << endl;

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

    cout << endl;

    return 0;
}
Wyniki
Równania czwartego stopnia typu
  4    3    2
ax + bx + cx + dx = 0
-------------------------------

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

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

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

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

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

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

Algorytm nr 3

Jeśli uda ci się znaleźć jeden pierwiastek x1 równania ax4+bx3+cx2+dx+e = 0, to wielomian tego równania możesz podzielić za pomocą schematu Hornera przez dwumian (x - x1):

W rezultacie otrzymasz równanie sześcienne, które rozwiązujesz metodami opisanymi w poprzednim rozdziale. Problem sprowadza się zatem do znalezienia pierwiastka rzeczywistego. Jeśli taki pierwiastek istnieje, to znajdziemy go jedną z metod przybliżonych, np. metodą Newtona. Przyjrzyjmy się wykresowi funkcji wielomianowej czwartego stopnia:

obrazek

W zależności od znaku współczynnika a mamy dwa wykresy, tutaj w kolorze niebieskim i czerwonym. Pierwsza pochodna tej funkcji wyraża się wzorem:

Ponieważ pochodna jest wielomianem sześciennym, to zawsze istnieje takie x0, dla którego pochodna f '( x0 ) = 0. Wszystkie takie punkty są punktami krytycznymi funkcji wyjściowej. Jeśli pochodna w danym punkcie jest zerowa, to styczna w tym punkcie do wykresu funkcji jest równoległa do osi OX i metoda Newtona załamuje się. Zatem poszukując pierwiastka musimy tak dobrać punkt startowy, aby na drodze do pierwiastka pierwsza pochodna nie zerowała się. Musimy więc najpierw znaleźć wszystkie punkty krytyczne i dokonać odpowiednich wyborów. Punkty krytyczne znajdziemy rozwiązując równanie:

Równanie to rozwiązujemy metodą z poprzedniego rozdziału. Jako rozwiązania interesują nas tylko pierwiastki rzeczywiste, których może być jeden lub trzy. Rozważmy te przypadki:

Jeden punkt krytyczny K

Funkcja wyjściowa ma w tym przypadku jedno maksimum lub jedno minimum w punkcie xK.

obrazek

Zwrot ramion wykresu zależy od znaku współczynnika a. Jeśli a jest dodatnie, to ramiona są skierowane do góry, a punkt K jest minimum funkcji (jak na obrazku powyżej). Jeśli a jest ujemne, to ramiona są skierowane w dół i punkt K jest maksimum funkcji:

obrazek

Sprawdzamy teraz wartość f ( xK ):

W takim przypadku xK jest pierwiastkiem równania wyjściowego i możemy pominąć metodę Newtona, przyjmując x1 = xK.

Mamy maksimum/minimum w xK i dwa pierwiastki po obu stronach xK. Za punkt startowy wybieramy x0 = xK – 1 lub x0 = xK + 1 i szukamy pierwiastka metodą Newtona.

W pozostałych przypadkach funkcja  f ( x ) nie ma pierwiastków rzeczywistych. Obliczenia przerywamy.

Trzy punkty krytyczne K1, K2 i K3

Funkcja wyjściowa  f ( x ) posiada minima/maksima lokalne:

obrazek obrazek

Porządkujemy punkty w kolejności:

Sprawdzamy, czy w którymś z tych punktów funkcja  f ( x ) przyjmuje wartość 0. Jeśli tak, to mamy jeden pierwiastek x1 i pomijamy metodę Newtona.

Sprawdzamy, czy w dwóch sąsiednich punktach xK1 i xK2 lub xK2 i xK3 wartości funkcji mają różne znaki. Jeśli tak, to za punkt startowy x0 wybieramy średnią arytmetyczną tych punktów i uruchamiamy metodę Newtona, aby znaleźć pierwiastek x1.

Sprawdzamy przypadki:

Za punkt startowy wybieramy x0 = xK1 1 i uruchamiamy metodę Newtona.

Za punkt startowy wybieramy x0 = xK3 + 1 i uruchamiamy metodę Newtona.

W przeciwnym razie funkcja  f ( x ) nie ma pierwiastków rzeczywistych i obliczenia przerywamy.

Metoda Newtona

Gdy zostanie ustalony punkt startowy x0, wykorzystujemy metodę Newtona do znalezienia pierwiastka rzeczywistego x1 funkcji wyjściowej  f ( x ).

Schemat Hornera

Gdy znany jest pierwiastek x1, dzielimy wielomian funkcji wyjściowej  f ( x ) przez dwumian (x - x1). Zadanie to wykonujemy przy pomocy schematu Hornera:

Równanie sześcienne

W rezultacie otrzymujemy wielomian sześcienny. Rozwiązujemy równanie a'x3+b'x2+c'x+d' = 0. Otrzymujemy pozostałe 3 pierwiastki x2, x3 i x4 i mamy pełne rozwiązanie równania czwartego stopnia.

Powyższa metoda pozwala znaleźć wszystkie pierwiastki równania czwartego stopnia, jeśli wśród nich występują pierwiastki rzeczywiste.

Metoda postępowania jest zatem następująca:

  1. Wyznaczamy punkty krytyczne funkcji  f ( x ) = ax4 + bx3 + cx2 + dx +e.
  2. Na podstawie punktów krytycznych sprawdzamy, czy wielomian funkcji ma pierwiastek rzeczywisty. Jeśli tak, to wyznaczamy go metodą Newtona. Jeśli nie, obliczenia przerywamy.
  3. Dzielimy za pomocą schematu Hornera wielomian funkcji przez dwumian x - x1, gdzie x1 jest pierwiastkiem wyznaczonym w punkcie 2. W wyniku otrzymujemy wielomian sześcienny.
  4. Rozwiązujemy równanie sześcienne i otrzymujemy pozostałe trzy pierwiastki.

Jak widzisz, metoda ta łączy ze sobą kilka opisanych wcześniej metod. Rozwiązanie otrzymamy, jeśli równanie wyjściowe posiada pierwiastek rzeczywisty. W przypadku samych pierwiastków zespolonych algorytm zawodzi. Należy wtedy zastosować jedną z metod ogólnych, które opisujemy w dalszej części tego rozdziału.

Przykładowy program wykorzystujący opisaną tutaj metodę jest następujący:

Przykładowy program w języku C++
// Równanie czwartego stopnia
// (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 oblicza wartość wielomianu metodą Hornera
//--------------------------------------------------

double f(double a, double b, double c, double d, double e, double x)
{
    double w = a;

    w = w * x + b;
    w = w * x + c;
    w = w * x + d;
    w = w * x + e;

    return w;
}
// Funkcja rozwiązująca równanie sześcienne
// a,b,c,d - współczynniki wielomianu
// cmplx - true, x1 = pierwiastek, x2 = re, x3 = +/- im
//         false, x1,x2,x3 = pierwiastki rzeczywiste
// x1,x2,x3 - pierwiastki
//-----------------------------------------
void ce(double a, double b, double c, double d,
        bool & cmplx,
        double & x1, double & x2, double & x3)
{
    double p,q,delta;

    cmplx = false;     // Znacznik wyniku zespolonego

    // 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;
    }
}

// Funkcja rozwiązująca równanie czwartego stopnia
// a,b,c,d,e - współczynniki wielomianu
//-----------------------------------------
void qe(double a, double b, double c, double d, double e)
{
    double x0,x1,x2,x3,x4;
    bool cmplx, noresult, newton;

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

    // Wyznaczamy pierwiastki pochodnej

    ce(4*a, 3*b, 2*c, d, cmplx, x1, x2, x3);
    noresult = false; // Znacznik braku pierwiastków
    newton = true;    // Znacznik obliczania pierwiastka metodą Newtona

    if(cmplx) // Jeden pierwiastek rzeczywisty
    {
        double yk = f(a,b,c,d,e,x1);        // Obliczamy wartość w punkcie xk
        if(fabs(yk) < epsy) newton = false; // Mamy jeden pierwiastek
        else if(yk * a < 0) x0 = x1 - 1;    // Punkt startowy dla metody Newtona
        else noresult = true;               // Brak pierwiastków rzeczywistych
    }
    else // Trzy pierwiastki rzeczywiste
    {
        // Porządkujemy punkty krytyczne

        if(x1 > x2) swap(x1,x2);
        if(x2 > x3) swap(x2,x3);
        if(x1 > x2) swap(x1,x2);

        double yk1,yk2,yk3;

        // Obliczamy wartości wielomianu w punktach krytycznych

        yk1 = f(a,b,c,d,e,x1);
        yk2 = f(a,b,c,d,e,x2);
        yk3 = f(a,b,c,d,e,x3);

        // Sprawdzamy, czy w którymś z punktów krytycznych wielomian zeruje się

        if(fabs(yk1) < epsy) newton = false;
        else if(fabs(yk2) < epsy)
        {
            newton = false;
            x1 = x2;
        }
        else if(fabs(yk3) < epsy)
        {
            newton = false;
            x1 = x3;
        }
        else if(yk1 * yk2 < 0) x0 = (x1 + x2) / 2; // Punkty startowe dla metody Newtona
        else if(yk2 * yk3 < 0) x0 = (x2 + x3) / 2;
        else if(yk1 * a < 0)   x0 = x1 - 1;
        else if(yk3 * a < 0)   x0 = x3 + 1;
        else noresult = true;
    }

    // Jeśli brak pierwiastków, przerywamy obliczenia

    if(noresult)
    {
        cout << "BRAK PIERWIASTKÓW RZECZYWISTYCH" << endl
             << "Obliczenia zostały przerwane!!!" << endl << endl;
        return;
    }

    // Jeśli nie znaleziono pierwiastka, uruchamiamy metodę Newtona

    if(newton)
    {
        double epsx = 1e-14; // Dokładność wyznaczania zera
        int n       = 64;    // Maksymalna liczba obiegów
        double f0,f1;
        bool result = false;

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

           f0 = f(a,b,c,d,e,x0);

           // Sprawdzamy, czy funkcja jest dostatecznie bliska zeru

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

           // Obliczamy wartość pierwszej pochodnej funkcji

           f1 = f(0,4*a,3*b,2*c,d,x0);

           // 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)
        {
            cout << "BŁĄD W METODZIE NEWTONA" << endl
                 << "Obliczenia zostały przerwane!!!" << endl << endl;
            return;
        }

        x1 = x0; // Pierwszy pierwiastek z metody Newtona
    }


    // W x1 mamy pierwiastek rzeczywisty
    // Schematem Hornera dzielimy wielomian równania przez dwumian (x - x1)

    double aa,bb,cc,dd;

    aa = a;
    bb = b + aa * x1;
    cc = c + bb * x1;
    dd = d + cc * x1;

    // Obliczamy pierwiastki wynikowego równania sześciennego

    ce(aa,bb,cc,dd,cmplx,x2,x3,x4);

    // Wyświetlamy wyniki

    cout << "x1 = " << x1 << endl
         << "x2 = " << x2 << endl;
    if(cmplx) cout << "x3 = " << x3 << " + " << -x4 << "i" << endl
                   << "x4 = " << x3 << " + " << x4 << "i";
    else      cout << "x3 = " << x3 << endl
                   << "x4 = " << x4;
    cout << endl;
}

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

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

    cout << fixed;
    cout << "Równania czwartego stopnia typu" << endl
         << "  4    3    2" << endl
         << "ax + bx + cx + dx + e = 0" << endl
         << "-------------------------------" << endl;

    qe(1,-10,35,-50,24);
    qe(1,-4,3,4,-4);
    qe(1,1,-7,-1,6);

    cout << endl;

    return 0;
}
Wynik
Równania czwartego stopnia typu
  4    3    2
ax + bx + cx + dx + e = 0
-------------------------------

Równanie: 1.000000 * x^4 + -10.000000 * x^3 + 35.000000 * x^2 + -50.000000 * x + 24.000000
x1 = 2.000000
x2 = 3.000000
x3 = 1.000000
x4 = 4.000000

Równanie: 1.000000 * x^4 + -4.000000 * x^3 + 3.000000 * x^2 + 4.000000 * x + -4.000000
x1 = 2.000000
x2 = 1.000000
x3 = -1.000000
x4 = 2.000000

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

Algorytm nr 4

Równanie czwartego stopnia można również rozwiązać algebraicznie za pomocą jednej z kilku metod. Tutaj przedstawimy metodę, która wykorzystuje pośrednie równanie sześcienne. Metoda pochodzi od włoskiego matematyka Lodovico Ferrari. Najpierw zobaczmy, jaką naturę mogą mieć pierwiastki równania czwartego stopnia:

Mamy ogólne równanie czwartego stopnia:

Załóżmy, że równanie to posiada pierwiastek rzeczywisty x1. W takim przypadku wielomian równania czwartego stopnia możemy zawsze podzielić przez dwumian (x - x1) :

Pozostałe 3 pierwiastki są pierwiastkami równania:

Ponieważ jest to równanie sześcienne, to pozostałe 3 pierwiastki mogą być albo wszystkie rzeczywiste, albo jeden rzeczywisty i dwa zespolone. Wynika z tego, iż równanie czwartego stopnia może mieć:

  • cztery pierwiastki rzeczywiste,
  • dwa pierwiastki rzeczywiste i dwa zespolone,
  • cztery pierwiastki zespolone.

Metoda rozwiązywania równania czwartego stopnia jest następująca:

Mamy ogólne równanie czwartego stopnia:

obrazek

Wyliczamy 3 kolejne współczynniki wg wzorów:

Z wyliczonych współczynników tworzymy pomocnicze równanie sześcienne:

Równanie to rozwiązujemy za pomocą metod opisanych w poprzednim rozdziale. Otrzymujemy trzy pierwiastki y1, y2 i y3. Wybieramy z nich dwa pierwiastki, które są różne od zera (jeśli mamy jeden pierwiastek rzeczywisty i dwa zespolone, zawsze wybieramy pierwiastki zespolone, które są liczbami sprzężonymi) i wyliczamy ich pierwiastki (mogą być zespolone):

Dodatkowo wyliczamy:

Pierwiastki równania wyjściowego otrzymujemy ze wzorów:

W działaniach mogą pojawiać się liczby zespolone, dlatego rachunki najlepiej przeprowadzać na liczbach zespolonych, dzięki czemu unikniemy rozważania różnych przypadków. Podobnie postąpiliśmy w algorytmie nr 1. Zasady działań na liczbach zespolonych są następujące (więcej na temat liczb zespolonych znajdziesz tutaj):

Poniższy program jest przykładową realizacją tego algorytmu:

Przykładowy program w języku C++
// Równanie czwartego stopnia
// (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;

// Liczby zespolone
//-----------------

// Definicja typu zespolonego
struct cmplx
{
    double re, im;
};

// Mnożenie liczb zespolonych
// c = a * b
//----------------------------
void cmplx_mul(cmplx a, cmplx b, cmplx & c)
{
    c.re = a.re * b.re - a.im * b.im;
    c.im = a.re * b.im + a.im * b.re;
}

// Dzielenie liczb zespolonych
// c = a / b
//----------------------------
void cmplx_div(cmplx a, cmplx b, cmplx & c)
{
    double d = b.re * b.re + b.im * b.im;
    c.re = (a.re * b.re + a.im * b.im) / d;
    c.im = (a.im * b.re - a.re * b.im) / d ;
}

// Pierwiastkowanie liczb zespolonych
// b = sqrt(a)
//--------------------------------------
void cmplx_sqrt(cmplx a, cmplx & b)
{
    double x,y;
    x = a.re;
    y = a.im;
    double c = sqrt(x * x + y * y);
    b.re = sqrt((c + x) / 2);
    b.im = sqrt((c - x) / 2);
    if(a.im < 0) b.im = - b.im;
}

// Wyświetlanie liczby zespolonej
void cmplx_print(cmplx a)
{
    cout << a.re;
    if(fabs(a.im) >= epsy)
    {
        if(a.im < 0) cout << " - ";
        else         cout << " + ";
        cout << fabs(a.im) << "i";
    }
    cout << endl;
}

// Funkcja rozwiązująca równanie sześcienne
// a,b,c,d - współczynniki wielomianu
// cmplx - true, x1 = pierwiastek, x2 = re, x3 = +/- im
//         false, x1,x2,x3 = pierwiastki rzeczywiste
// x1,x2,x3 - pierwiastki
//-----------------------------------------
void ce(double a, double b, double c, double d,
        bool & cmpl,
        double & x1, double & x2, double & x3)
{
    double p,q,delta;

    cmpl = false;     // Znacznik wyniku zespolonego

    // 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
    {
        cmpl = 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;
    }
}

// Funkcja rozwiązująca równanie czwartego stopnia
// a,b,c,d,e - współczynniki wielomianu
//-----------------------------------------
void qe(double a, double b, double c, double d, double e)
{
    cout << endl
         << "Równanie: "
         << a << " * x^4 + "
         << b << " * x^3 + "
         << c << " * x^2 + "
         << d << " * x + "
         << e << endl;
    double f,g,h;

    // Wyliczamy współczynniki

    f = (c - 3 * b * b / 8 / a) / a;
    g = (d + b * b * b / 8 / a / a - b * c / 2 / a) / a;
    h = (e - 3 * b * b * b * b / 256 / a / a / a + b * b * c / 16 / a / a - b * d / 4 / a) / a;

    double x1,x2,x3;
    bool cmpl;

    // Tworzymy pomocnicze równanie sześcienne i wyliczamy jego pierwiastki

    ce(1,f / 2, (f * f - 4 * h) / 16, - g * g / 64, cmpl, x1, x2, x3);

    cmplx y1, y2;

    if(cmpl) // Pierwiastki zespolone
    {
        y1.re = y2.re = x2;
        y1.im = -x3;
        y2.im = x3;
    }
    else  // Wybieramy dwa niezerowe pierwiastki
    {
        y1.im = y2.im = 0;
        if(fabs(x1) < epsy)
        {
            y1.re = x2;
            y2.re = x3;
        }
        else if(fabs(x2) < epsy)
        {
            y1.re = x1;
            y2.re = x3;
        }
        else
        {
            y1.re = x1;
            y2.re = x2;
        }
    }

    // Obliczamy ich pierwiastki zespolone

    cmplx p,q,r,t,u;

    double s;

    cmplx_sqrt(y1,p);
    cmplx_sqrt(y2,q);

    // Wyliczamy współczynniki s i r

    s = b / 4 / a;

    cmplx_mul(p,q,t);
    u.re = - g / 8;
    u.im = 0;
    cmplx_div(u,t,r);

    // Obliczamy pierwiastki zespolone równania czwartego stopnia

    cmplx r1,r2,r3,r4;

    r1.re = p.re + q.re + r.re - s;
    r1.im = p.im + q.im + r.im;

    r2.re = p.re - q.re - r.re - s;
    r2.im = p.im - q.im - r.im;

    r3.re = -p.re + q.re - r.re - s;
    r3.im = -p.im + q.im - r.im;

    r4.re = -p.re - q.re + r.re - s;
    r4.im = -p.im - q.im + r.im;

    // Wyświetlamy wyniki

    cout << "x1 = "; cmplx_print(r1);
    cout << "x2 = "; cmplx_print(r2);
    cout << "x3 = "; cmplx_print(r3);
    cout << "x4 = "; cmplx_print(r4);
}

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

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

    cout << fixed;
    cout << "Równania czwartego stopnia typu" << endl
         << "  4    3    2" << endl
         << "ax + bx + cx + dx + e = 0" << endl
         << "-------------------------------" << endl;

    qe(1,-10,35,-50,24);
    qe(1,-4,3,4,-4);
    qe(1,1,-7,-1,6);
    qe(5,5,5,5,5);

    cout << endl;

    return 0;
}
Wynik
Równania czwartego stopnia typu
  4    3    2
ax + bx + cx + dx + e = 0
-------------------------------

Równanie: 1.000000 * x^4 + -10.000000 * x^3 + 35.000000 * x^2 + -50.000000 * x + 24.000000
x1 = 4.000000
x2 = 2.000000
x3 = 3.000000
x4 = 1.000000

Równanie: 1.000000 * x^4 + -4.000000 * x^3 + 3.000000 * x^2 + 4.000000 * x + -4.000000
x1 = 2.000000
x2 = 2.000000
x3 = 1.000000
x4 = -1.000000

Równanie: 1.000000 * x^4 + 1.000000 * x^3 + -7.000000 * x^2 + -1.000000 * x + 6.000000
x1 = -1.000000
x2 = 2.000000
x3 = 1.000000
x4 = -3.000000

Równanie: 5.000000 * x^4 + 5.000000 * x^3 + 5.000000 * x^2 + 5.000000 * x + 5.000000
x1 = 0.309017 + 0.951057i
x2 = -0.809017 - 0.587785i
x3 = -0.809017 + 0.587785i
x4 = 0.309017 - 0.951057i
Na początek:  podrozdziału   strony 

Równania wyższych stopni

Równanie stopnia czwartego jest ostatnim równaniem wielomianowym, dla którego istnieje algebraiczna metoda rozwiązywania. Udowodniono, że nie ma takiej metody dla równań o stopniu wyższym niż 4. Nie oznacza to oczywiście, że równania takie nie posiadają pierwiastków. Nie można jedynie znaleźć tych pierwiastków ogólną metodą o skończonej liczbie operacji. Pozostają nam zatem metody przybliżone lub metody rozwiązujące szczególne przypadki tych równań.

Do znalezienia rozwiązania z dowolną dokładnością można wykorzystać metodę Newtona oraz podział wielomianu schematem Hornera, gdy uda nam się znaleźć pierwiastek przybliżony. W efekcie otrzymamy równanie niższego stopnia. Korzystamy tutaj z tej samej procedury, którą stosowaliśmy dla równań 3-go i 4-go stopnia:

  1. Wyznaczamy punkty krytyczne.
  2. Porządkujemy wyznaczone punkty krytyczne rosnąco.
  3. Sprawdzamy, czy pierwiastek leży w jednym z punktów krytycznych. Jeśli tak, to przechodzimy do kroku 6.
  4. Sprawdzamy, czy w dwóch sąsiednich punktach krytycznych funkcja ma różne znaki. Jeśli tak, to pierwiastek znajduje się pomiędzy tymi punktami krytycznymi, znajdujemy go metodą Newtona i przechodzimy do kroku 6.
  5. Sprawdzamy, czy pierwiastek może być po lewej stronie pierwszego punktu krytycznego lub po prawej stronie ostatniego punktu krytycznego. Jeśli tak, to znajdujemy go metodą Newtona i przechodzimy do kroku 6. Jeśli nie, to równanie nie ma pierwiastków rzeczywistych i przerywamy obliczenia.
  6. Mamy wyznaczony pierwszy pierwiastek x1. Wielomian równania dzielimy za pomocą schematu Hornera przez dwumian (x - x1) i otrzymujemy równanie wynikowe o stopniu niższym o 1.
  7. Rozwiązujemy równanie wynikowe i otrzymujemy pełne rozwiązanie równania wyjściowego.

Zwróć uwagę, iż takie podejście pozwala rozwiązać równanie stopnia n, jeśli potrafimy rozwiązać równanie stopnia n –1. Dodatkowo pamiętaj, iż równania o stopniach nieparzystych posiadają zawsze co najmniej jeden pierwiastek rzeczywisty.

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.