Zagadnienie przydziału i algorytm węgierski — model, zapis liniowy i rozwiązanie krok po kroku.

Mgr inż.

Sebastian Dziarmaga-Działyński

Wirtualna kancelaria korepetytorska i konsultacyjna – wszechwiedza.pl

Zagadnienie przydziału polega na optymalnym przypisaniu wykonawców do zadań. Każdy wykonawca ma otrzymać jedno zadanie, a każde zadanie ma zostać wykonane przez dokładnie jednego wykonawcę. Celem może być minimalizacja łącznego kosztu lub czasu albo maksymalizacja łącznego zysku, efektu czy korzyści. Jedną z najpopularniejszych metod rozwiązywania takich problemów minimalizacyjnych jest algorytm węgierski.

Na czym polega zagadnienie przydziału?

Zagadnienie przydziału jest jednym z podstawowych modeli badań operacyjnych. W najprostszym wariancie mamy taką samą liczbę wykonawców i zadań. Każdemu wykonawcy można przypisać każde zadanie, ale koszt, czas lub korzyść zależą od konkretnej pary.

Przykładowo można przydzielać:

• pracowników do stanowisk pracy,
• kierowców do tras,
• maszyny do zleceń produkcyjnych,
• techników do zgłoszeń serwisowych,
• wykładowców do grup zajęciowych,
• programistów do modułów projektu,
• pojazdy do rejonów dostaw.

W każdej z tych sytuacji obowiązuje zasadnicza reguła: jeden wykonawca otrzymuje jedno zadanie, a jedno zadanie zostaje powierzone jednemu wykonawcy. Nie chodzi więc o podział pracy w dowolnych proporcjach, lecz o wybór najlepszego pełnego skojarzenia dwóch zbiorów.

Macierz kosztów i zadanie niezrównoważone

Dane do zadania przydziału zapisuje się zwykle w postaci macierzy. Wiersze reprezentują wykonawców, kolumny reprezentują zadania, a liczba w komórce oznacza koszt wykonania określonego zadania przez określonego wykonawcę.

Symbol \(c_{ij}\) oznacza koszt przydzielenia wykonawcy \(i\) do zadania \(j\). Zamiast kosztu można w ten sam sposób zapisać czas realizacji, zużycie materiału, emisję, liczbę przejechanych kilometrów albo inną wielkość, którą chcemy minimalizować.

Zadanie zrównoważone

Klasyczna postać algorytmu węgierskiego dotyczy macierzy kwadratowej, czyli sytuacji, w której liczba wykonawców jest równa liczbie zadań. Jeśli mamy czterech wykonawców i cztery zadania, otrzymujemy macierz \(4\times4\).

Zadanie niezrównoważone

W praktyce liczba wykonawców i zadań często nie jest taka sama. Gdy wykonawców jest mniej niż zadań, dodajemy fikcyjnych wykonawców. Gdy wykonawców jest więcej, dodajemy fikcyjne zadania. W ten sposób uzupełniamy macierz do postaci kwadratowej.

Jeżeli brak realizacji zadania nie powoduje dodatkowego kosztu, w komórkach związanych z fikcyjnym wykonawcą lub fikcyjnym zadaniem można wpisać zero. Gdy jednak pominięcie zadania oznacza stratę, karę umowną lub koszt alternatywny, należy wpisać wartość odzwierciedlającą tę konsekwencję.

Model programowania liniowego ze zmiennymi binarnymi

Zagadnienie przydziału można zapisać jako model programowania liniowego. Wprowadzamy zmienne decyzyjne:

Jeżeli minimalizujemy łączny koszt, funkcja celu ma postać:

Każdy wykonawca ma otrzymać dokładnie jedno zadanie:

Każde zadanie ma zostać przydzielone dokładnie jednemu wykonawcy:

Ostatecznie określamy dziedzinę zmiennych:

Model można rozwiązać metodą programowania całkowitoliczbowego, na przykład w Excel Solverze lub w specjalistycznym oprogramowaniu optymalizacyjnym. Dla klasycznego problemu przydziału bardzo wygodny jest jednak algorytm węgierski, ponieważ wykorzystuje specyficzną strukturę tego zadania.

Minimalizacja a maksymalizacja

Algorytm węgierski najczęściej przedstawia się dla zagadnień minimalizacyjnych. Jeżeli w macierzy znajdują się zyski, efekty lub korzyści, które chcemy maksymalizować, można najpierw przekształcić zadanie maksymalizacyjne do równoważnego zadania minimalizacyjnego.

Niech \(p_{ij}\) oznacza korzyść uzyskaną przy przypisaniu wykonawcy \(i\) do zadania \(j\). Najpierw wyznaczamy największy element całej macierzy:

Następnie dla każdej komórki tworzymy koszt przekształcony:

Największe korzyści zmieniają się wtedy w najmniejsze koszty. Minimalizacja sumy wartości \(c_{ij}\) prowadzi do tego samego przydziału, który wcześniej maksymalizował sumę korzyści \(p_{ij}\).

Algorytm węgierski

Algorytm węgierski prowadzi do rozwiązania przez kolejne przekształcanie macierzy kosztów. Kluczowe jest to, że po redukcjach nadal otrzymujemy przydział optymalny dla pierwotnego problemu.

Krok 1. Redukcja wierszy

W każdym wierszu odnajdujemy najmniejszy element i odejmujemy go od wszystkich wartości w tym wierszu. W każdym wierszu pojawi się wtedy przynajmniej jedno zero.

Krok 2. Redukcja kolumn

W każdej kolumnie z macierzy otrzymanej po redukcji wierszy odnajdujemy najmniejszy element i odejmujemy go od wszystkich wartości w tej kolumnie. Po tej operacji w każdej kolumnie występuje co najmniej jedno zero.

Krok 3. Pokrycie zer minimalną liczbą linii

Wszystkie zera w macierzy pokrywamy minimalną liczbą linii poziomych i pionowych. Jeżeli liczba użytych linii jest równa wymiarowi macierzy, można przejść do wyboru niezależnych zer. Jeżeli liczba linii jest mniejsza, potrzebna jest korekta macierzy.

Krok 4. Korekta macierzy

Gdy liczba linii pokrywających zera jest mniejsza niż liczba wierszy, znajdujemy najmniejszy element niepokryty liniami. Następnie:

1. odejmujemy tę wartość od wszystkich elementów niepokrytych,
2. dodajemy ją do elementów leżących na przecięciu dwóch linii,
3. nie zmieniamy elementów pokrytych dokładnie jedną linią.

Po korekcie ponownie pokrywamy wszystkie zera minimalną liczbą linii. Czynność powtarzamy aż do momentu, gdy liczba linii będzie równa wymiarowi macierzy.

Krok 5. Wybór niezależnych zer

Niezależne zera to takie zera, z których żadne dwa nie leżą w tym samym wierszu ani w tej samej kolumnie. Wybieramy po jednym zerze dla każdego wykonawcy i dla każdego zadania. Wybrane zera określają ostateczny przydział.

Przykład rozwiązany krok po kroku

Załóżmy, że czterech pracowników \(A\), \(B\), \(C\) i \(D\) ma zostać przydzielonych do czterech zadań \(Z_1\), \(Z_2\), \(Z_3\) i \(Z_4\). Liczby w tabeli oznaczają koszt wykonania danego zadania przez danego pracownika.

Krok 1. Redukcja wierszy

Od każdego elementu wiersza \(A\) odejmujemy 69, od każdego elementu wiersza \(B\) odejmujemy 37, od każdego elementu wiersza \(C\) odejmujemy 5, a od każdego elementu wiersza \(D\) odejmujemy 8.

Krok 2. Redukcja kolumn

Po redukcji wierszy minima kolumn wynoszą kolejno: \(0\), \(0\), \(0\) oraz \(15\). Tylko w czwartej kolumnie trzeba więc odjąć wartość 15 od wszystkich elementów.

Krok 3. Pokrycie zer liniami

Wszystkie zera można pokryć trzema liniami: linią poziomą wiersza \(D\) oraz liniami pionowymi kolumn \(Z_2\) i \(Z_3\). Ponieważ macierz ma wymiar \(4\times4\), a liczba linii jest równa 3, nie można jeszcze wybrać czterech niezależnych zer. Trzeba wykonać korektę macierzy.

Najmniejszym elementem niepokrytym liniami jest liczba \(6\). Odejmujemy ją od wszystkich elementów niepokrytych, dodajemy ją do elementów leżących na przecięciu dwóch linii, a pozostałe elementy pozostawiamy bez zmian.

Krok 4. Korekta macierzy

Możemy teraz wybrać cztery niezależne zera:

• wykonawca \(A\) otrzymuje zadanie \(Z_3\),
• wykonawca \(B\) otrzymuje zadanie \(Z_2\),
• wykonawca \(C\) otrzymuje zadanie \(Z_1\),
• wykonawca \(D\) otrzymuje zadanie \(Z_4\).

Obliczając wartość funkcji celu, wracamy do macierzy pierwotnej. Zera z macierzy zredukowanej służą wyłącznie do wskazania optymalnego układu przydziałów.

Najczęstsze błędy

• Odczytywanie kosztu z macierzy zredukowanej. Wartość funkcji celu liczymy zawsze na podstawie pierwotnej macierzy kosztów.
• Wybieranie dwóch zer z jednego wiersza lub jednej kolumny. Każdy wykonawca i każde zadanie muszą wystąpić dokładnie raz.
• Pomijanie elementu fikcyjnego w zadaniu niezrównoważonym. Algorytm klasyczny wymaga macierzy kwadratowej.
• Uznawanie, że każde zero powinno zostać wybrane. Wybieramy tylko niezależne zera, czyli takie, które nie kolidują ze sobą w wierszach i kolumnach.
• Niepoprawna korekta macierzy. Najmniejszą wartość niepokrytą odejmujemy od elementów niepokrytych, ale dodajemy ją wyłącznie na przecięciach dwóch linii.
• Niewłaściwe przekształcenie maksymalizacji. Przed uruchomieniem algorytmu trzeba zamienić macierz korzyści na macierz kosztów przekształconych.

Podsumowanie

Zagadnienie przydziału służy do wyboru najlepszego pełnego przypisania wykonawców do zadań. Można je formalnie zapisać jako problem programowania liniowego ze zmiennymi binarnymi, a w klasycznej wersji minimalizacyjnej wygodnie rozwiązać algorytmem węgierskim.

Najważniejsze etapy algorytmu to redukcja wierszy, redukcja kolumn, pokrywanie zer minimalną liczbą linii, ewentualna korekta macierzy oraz wybór niezależnych zer. Warto pamiętać, że zredukowana macierz pomaga odnaleźć optymalny przydział, ale całkowity koszt lub zysk obliczamy zawsze na podstawie danych początkowych.

Utworzono: 20.06.2026 | Zmodyfikowano: 21.06.2026

Kinematyka ciała sztywnego – ruch postępowy, obrotowy i płaski

Kinematyka ciała sztywnego opisuje zmianę położenia i orientacji obiektu, którego odkształcenia można pominąć. W odróżnieniu od punktu materialnego nie wystarczy określić ruchu jednego punktu — trzeba również uwzględnić obrót całego ciała oraz zależności między prędkościami i przyspieszeniami jego poszczególnych punktów.

Do podstawowych rodzajów ruchu ciała sztywnego należą ruch postępowy, ruch obrotowy wokół stałej osi oraz ruch płaski, będący złożeniem przesunięcia i obrotu. Ważnymi przykładami ruchu płaskiego są toczenie koła bez poślizgu oraz ruch członów mechanizmów.

Opis ruchu ciała sztywnego wykorzystuje wiadomości przedstawione wcześniej w artykułach Kinematyka punktu materialnego – równania ruchu, prędkość, przyspieszenie i promień krzywizny oraz Ruch po okręgu – prędkość kątowa, przyspieszenie dośrodkowe, okres i częstotliwość.

Czym jest ciało sztywne?

Ciało sztywne jest idealizowanym układem punktów materialnych, w którym odległość między dowolnymi dwoma punktami nie zmienia się podczas ruchu.

Dla dowolnych punktów \(A\) i \(B\) należących do ciała warunek sztywności ma postać:

Równoważnie można zapisać:

gdzie:

jest wektorem położenia punktu \(B\) względem punktu \(A\).

Warunek ten oznacza, że ciało może zmieniać położenie i orientację, ale w przyjętym modelu nie zmienia kształtu ani wymiarów.

Punkt materialny a ciało sztywne

W modelu punktu materialnego pomijamy rozmiary, kształt i orientację obiektu. Jego ruch opisujemy za pomocą jednego wektora położenia:

W przypadku ciała sztywnego trzeba określić nie tylko położenie wybranego punktu odniesienia, lecz także orientację całego ciała.

W ruchu płaskim opis może obejmować:

gdzie \(A\) jest wybranym punktem ciała, a \(\varphi\) określa jego orientację w płaszczyźnie.

Położenie i orientacja ciała

Położenie ciała sztywnego można opisać przez wskazanie położenia jednego punktu odniesienia oraz orientacji układu współrzędnych związanego z ciałem.

Dla dwóch punktów \(A\) i \(B\):

Wektor \(\vec r_{B/A}\) ma stałą długość, ale podczas obrotu może zmieniać kierunek względem nieruchomego układu odniesienia.

Zmiana położenia punktu \(A\) odpowiada za przesunięcie ciała, natomiast zmiana kierunku wektora \(\vec r_{B/A}\) odpowiada za obrót.

Stopnie swobody

Liczba stopni swobody określa liczbę niezależnych współrzędnych potrzebnych do jednoznacznego opisania położenia ciała.

Ciało sztywne w przestrzeni

Swobodne ciało sztywne w przestrzeni ma sześć stopni swobody:

• trzy przesunięcia w kierunkach osi \(x\), \(y\) i \(z\);
• trzy niezależne obroty określające orientację ciała.

Ciało wykonujące ruch płaski

W ruchu płaskim ciało ma trzy stopnie swobody:

Są to dwa przesunięcia w płaszczyźnie oraz jeden obrót wokół osi prostopadłej do tej płaszczyzny.

Podstawowe rodzaje ruchu ciała sztywnego

W kinematyce wyróżniamy przede wszystkim:

• ruch postępowy — orientacja ciała nie zmienia się;
• ruch obrotowy wokół stałej osi — punkty ciała poruszają się po okręgach wokół osi;
• ruch płaski — ciało jednocześnie przesuwa się i obraca.

Ruch postępowy

W ruchu postępowym każdy odcinek związany z ciałem pozostaje równoległy do swojego wcześniejszego położenia. Orientacja ciała nie zmienia się.

Wektor położenia punktu \(B\) możemy zapisać jako:

W ruchu postępowym wektor \(\vec r_{B/A}\) jest stały nie tylko co do długości, lecz także co do kierunku. Po zróżniczkowaniu względem czasu otrzymujemy:

Po ponownym zróżniczkowaniu:

Wszystkie punkty ciała wykonującego ruch postępowy mają więc w tej samej chwili jednakowe prędkości i przyspieszenia.

Przykładami mogą być:

• kabina windy poruszająca się po pionowej prowadnicy;
• tłok przesuwający się w cylindrze;
• platforma pozostająca poziomo podczas ruchu po zakrzywionym torze;
• szuflada wysuwana z mebla.

Ruch obrotowy wokół stałej osi

W ruchu obrotowym wokół stałej osi wszystkie punkty ciała znajdujące się poza osią poruszają się po okręgach. Środki tych okręgów leżą na osi obrotu, a ich płaszczyzny są do niej prostopadłe.

Punkty należące do osi obrotu pozostają nieruchome:

Położenie kątowe ciała opisuje funkcja:

Prędkość kątowa wynosi:

Przyspieszenie kątowe:

Wszystkie punkty ciała mają w danej chwili tę samą prędkość kątową \(\vec\omega\) i to samo przyspieszenie kątowe \(\vec\varepsilon\).

Prędkość punktu w ruchu obrotowym

Dla punktu \(P\), którego położenie względem wybranego punktu osi opisuje wektor \(\vec r\), prędkość wynosi:

Wektor prędkości jest:

• styczny do okręgu opisywanego przez punkt;
• prostopadły do promienia prowadzącego od osi obrotu;
• zgodny ze zwrotem wyznaczonym przez iloczyn wektorowy.

Wartość prędkości liniowej punktu oddalonego od osi o odległość \(r\) wynosi:

Im dalej punkt znajduje się od osi, tym większa jest jego szybkość liniowa. Jeżeli:

to:

Przyspieszenie punktu w ruchu obrotowym

Przyspieszenie punktu obracającego się ciała składa się z dwóch wzajemnie prostopadłych części:

• przyspieszenia stycznego, związanego ze zmianą wartości prędkości kątowej;
• przyspieszenia normalnego, związanego ze zmianą kierunku prędkości.

Przyspieszenie styczne:

Jego wartość wynosi:

Przyspieszenie normalne:

Jest ono skierowane ku osi obrotu, a jego wartość wynosi:

Całkowite przyspieszenie punktu:

Ponieważ składowe styczna i normalna są prostopadłe, wartość przyspieszenia wynosi:

Ruch płaski ciała sztywnego

Ruch płaski zachodzi wtedy, gdy wszystkie punkty ciała poruszają się w płaszczyznach równoległych do jednej ustalonej płaszczyzny.

Ruch ten można przedstawić jako złożenie:

• ruchu postępowego wybranego punktu odniesienia;
• obrotu całego ciała wokół osi przechodzącej przez ten punkt i prostopadłej do płaszczyzny ruchu.

Wybór punktu odniesienia jest dowolny. Najczęściej przyjmuje się środek masy, punkt połączenia członów mechanizmu albo punkt, którego prędkość jest znana.

Zależność między prędkościami dwóch punktów

Dla dwóch punktów \(A\) i \(B\) należących do tego samego ciała sztywnego:

Prędkość punktu \(B\) jest sumą:

• prędkości wybranego punktu odniesienia \(A\);
• prędkości punktu \(B\) względem \(A\), wynikającej z obrotu ciała.

Składnik względny:

jest prostopadły do odcinka \(AB\), a jego wartość wynosi:

Zależność między przyspieszeniami dwóch punktów

Przyspieszenia dwóch punktów należących do ciała sztywnego są związane równaniem:

Składniki tego równania oznaczają:

• \(\vec a_A\) — przyspieszenie punktu odniesienia \(A\);
• \(\vec\varepsilon\times\vec r_{B/A}\) — względne przyspieszenie styczne punktu \(B\);
• \(\vec\omega\times(\vec\omega\times\vec r_{B/A})\) — względne przyspieszenie normalne, skierowane od \(B\) ku \(A\).

W ruchu płaskim wartości składowych względnych wynoszą:

Zależność przyspieszeń jest bardziej złożona niż zależność prędkości, ponieważ występują w niej dwa składniki związane z obrotem.

Chwilowy środek obrotu

W ruchu płaskim, jeżeli prędkość kątowa nie jest równa zeru, można zwykle wskazać punkt \(I\), którego prędkość w danej chwili jest równa zeru:

Punkt ten nazywamy chwilowym środkiem obrotu. W rozpatrywanej chwili ruch ciała można interpretować jako czysty obrót wokół punktu \(I\).

Dla dowolnego punktu \(A\):

Wartość prędkości:

Wektor prędkości punktu \(A\) jest prostopadły do odcinka \(IA\):

Analogicznie:

Jak wyznaczyć chwilowy środek obrotu?

Jeżeli znamy kierunki prędkości dwóch punktów ciała:

• przez pierwszy punkt prowadzimy prostą prostopadłą do jego prędkości;
• przez drugi punkt prowadzimy prostą prostopadłą do jego prędkości;
• punkt przecięcia tych prostych jest chwilowym środkiem obrotu.

Chwilowy środek obrotu może:

• leżeć wewnątrz ciała;
• leżeć na jego brzegu;
• znajdować się poza ciałem;
• zmieniać położenie wraz z czasem.

Toczenie koła bez poślizgu

Toczenie koła jest przykładem ruchu płaskiego. Koło jednocześnie:

• przesuwa się wraz ze swoim środkiem \(C\);
• obraca się wokół środka z prędkością kątową \(\omega\).

Jeżeli koło toczy się bez poślizgu po nieruchomym podłożu, punkt styku \(P\) ma chwilowo zerową prędkość:

Warunek toczenia bez poślizgu ma postać:

gdzie:

• \(v_C\) — szybkość środka koła;
• \(\omega\) — wartość prędkości kątowej;
• \(R\) — promień koła.

Dla koła poruszającego się w prawo obrót odbywa się zgodnie z ruchem wskazówek zegara.

Prędkości charakterystycznych punktów koła

Prędkość dowolnego punktu \(A\) koła jest sumą prędkości środka i prędkości wynikającej z obrotu:

Środek koła

Punkt styku z podłożem

Prędkość wynikająca z obrotu ma taką samą wartość jak prędkość środka, ale przeciwny zwrot:

Najwyższy punkt koła

Oba składniki prędkości mają zgodny zwrot:

Ponieważ:

otrzymujemy:

Podsumowując:

Przykłady obliczeniowe

Przykład 1. Prędkości punktów obracającej się tarczy

Tarcza obraca się z prędkością kątową:

Punkt \(A\) znajduje się w odległości:

a punkt \(B\) w odległości:

Szybkość punktu \(A\):

Szybkość punktu \(B\):

Oba punkty mają tę samą prędkość kątową, ale punkt \(B\), położony dalej od osi, ma większą szybkość liniową.

Przykład 2. Przyspieszenie punktu tarczy

Tarcza ma:

Rozpatrywany punkt znajduje się w odległości:

Przyspieszenie styczne:

Przyspieszenie normalne:

Wartość przyspieszenia całkowitego:

Przykład 3. Prędkość drugiego punktu ciała

Prędkość punktu \(A\) wynosi:

Prędkość kątowa ciała:

Położenie punktu \(B\) względem \(A\):

Składnik wynikający z obrotu:

Prędkość punktu \(B\):

Jej wartość wynosi:

Przykład 4. Koło toczące się bez poślizgu

Koło o promieniu:

toczy się bez poślizgu, a szybkość jego środka wynosi:

Prędkość kątowa:

Szybkość najwyższego punktu:

Szybkość punktu styku:

Jednostki i kontrola wymiarów

Kontrola wymiarów zależności \(v=\omega r\):

Radian jest w analizie wymiarowej traktowany jako jednostka bezwymiarowa.

Kontrola przyspieszenia stycznego:

Kontrola przyspieszenia normalnego:

Najczęstsze błędy

• traktowanie ciała sztywnego jak pojedynczego punktu materialnego;
• pomijanie zmiany orientacji ciała;
• uznawanie, że ruch postępowy zawsze musi być prostoliniowy;
• przyjmowanie jednakowych prędkości liniowych wszystkich punktów obracającego się ciała;
• mylenie wspólnej prędkości kątowej z prędkościami liniowymi punktów;
• rysowanie wektora prędkości punktu obracającego się ciała wzdłuż promienia zamiast stycznie do toru;
• pomijanie prostopadłości \(\vec v\perp\vec r\) w ruchu obrotowym;
• stosowanie \(v=\omega/r\) zamiast \(v=\omega r\);
• mylenie przyspieszenia stycznego z normalnym;
• kierowanie przyspieszenia normalnego na zewnątrz zamiast ku osi obrotu;
• dodawanie wartości prędkości punktów zamiast wykonywania dodawania wektorowego;
• pomijanie składnika normalnego w równaniu przyspieszeń dwóch punktów;
• uznawanie, że chwilowy środek obrotu musi znajdować się wewnątrz ciała;
• uznawanie, że chwilowy środek obrotu ma zawsze zerowe przyspieszenie;
• uznawanie punktu styku toczącego się koła za punkt trwale nieruchomy;
• przyjmowanie \(v_P=v_C\) dla punktu styku koła toczącego się bez poślizgu;
• przyjmowanie, że najwyższy punkt koła porusza się z szybkością \(v_C\) zamiast \(2v_C\);
• mylenie toczenia bez poślizgu z ruchem koła ślizgającego się po podłożu;
• niezgodność zwrotu obrotu koła ze zwrotem ruchu jego środka.

Podsumowanie

Warunek sztywności ciała ma postać:

W ruchu postępowym wszystkie punkty mają jednakowe prędkości i przyspieszenia:

W ruchu obrotowym wokół stałej osi:

W ruchu płaskim zależność między prędkościami dwóch punktów wynosi:

Dla przyspieszeń:

Jeżeli \(I\) jest chwilowym środkiem obrotu:

Podczas toczenia koła bez poślizgu:

Kinematyka ciała sztywnego pozwala opisywać ruch kół, tarcz, dźwigni, korb, tłoków oraz całych mechanizmów. Kluczowe znaczenie ma przy tym rozróżnienie ruchu całego ciała od ruchu jego poszczególnych punktów.

Numeryczne rozwiązywanie równań różniczkowych – metody Eulera i Rungego – Kutty

Mgr inż.

Sebastian Dziarmaga-Działyński

Wirtualna kancelaria korepetytorska i konsultacyjna – wszechwiedza.pl

Równania różniczkowe opisują zmiany zachodzące w czasie lub przestrzeni. Występują między innymi w mechanice, elektrotechnice, termodynamice, biologii, ekonomii, automatyce i teorii sterowania. Pozwalają modelować ruch ciał, przepływ ciepła, drgania, rozwój populacji, przebieg reakcji chemicznych oraz działanie obwodów elektrycznych.

Nie każde równanie różniczkowe można jednak rozwiązać analitycznie. Nawet jeżeli rozwiązanie dokładne istnieje, jego wyznaczenie może być bardzo trudne albo otrzymana postać może być mało użyteczna w praktycznych obliczeniach.

W takich sytuacjach stosuje się numeryczne metody rozwiązywania równań różniczkowych. Zamiast szukać jawnego wzoru funkcji, wyznaczamy jej przybliżone wartości w kolejnych punktach siatki.

W artykule omówimy między innymi:

• jawną i niejawną metodę Eulera,
• ulepszoną metodę Eulera, nazywaną metodą Heuna,
• metodę punktu środkowego,
• klasyczną metodę Rungego – Kutty czwartego rzędu,
• metody adaptacyjne,
• metody wielokrokowe,
• stabilność i równania sztywne.

Podstawowe znaczenie algorytmu numerycznego, zbieżności, stabilności i błędu zostało przedstawione w artykułach Metody numeryczne – czym są i do czego służą? oraz Błędy obliczeń numerycznych – rodzaje, źródła i ocena dokładności.

Zagadnienie początkowe

Rozważmy równanie różniczkowe pierwszego rzędu:

\[ y'(t)=f(t,y). \]

Samo równanie może mieć całą rodzinę rozwiązań. Aby wybrać jedną konkretną funkcję, podajemy warunek początkowy:

\[ y(t_0)=y_0. \]

Otrzymujemy w ten sposób zagadnienie początkowe, nazywane także zagadnieniem Cauchy’ego:

\[ \begin{cases} y'(t)=f(t,y),\\[2mm] y(t_0)=y_0. \end{cases} \]

Funkcja \(f(t,y)\) określa nachylenie rozwiązania w punkcie \((t,y)\). Jeżeli znamy aktualny stan układu, równanie różniczkowe mówi nam, w jakim kierunku rozwiązanie powinno się lokalnie zmieniać.

Rozwiązanie analityczne a rozwiązanie numeryczne

W rozwiązaniu analitycznym próbujemy wyznaczyć funkcję:

\[ y=y(t), \]

która spełnia równanie różniczkowe i warunek początkowy na całym rozpatrywanym przedziale.

Przykładowo dla zagadnienia:

\[ y'=-2y, \qquad y(0)=1 \]

rozwiązaniem dokładnym jest funkcja:

\[ y(t)=e^{-2t}. \]

Rozwiązanie numeryczne ma zwykle inną postać. Wybieramy punkty:

\[ t_0,\ t_1,\ t_2,\ldots,t_N \]

i obliczamy wartości:

\[ y_0,\ y_1,\ y_2,\ldots,y_N, \]

które powinny przybliżać dokładne rozwiązanie:

\[ y_n\approx y(t_n). \]

Podstawowym wynikiem obliczeń jest więc zbiór punktów:

\[ (t_0,y_0),\ (t_1,y_1),\ldots,(t_N,y_N). \]

Punkty można następnie połączyć odcinkami albo zastosować interpolację, aby otrzymać przybliżony wykres ciągły. Nadal jednak nie jest to na ogół symboliczny wzór funkcji, lecz jej numeryczna reprezentacja.

Siatka obliczeniowa i krok metody

W najprostszym przypadku stosujemy stały krok:

\[ h>0. \]

Kolejne punkty siatki wyznaczamy ze wzoru:

\[ t_{n+1}=t_n+h. \]

Stąd:

\[ t_n=t_0+nh. \]

Jeżeli chcemy obliczać rozwiązanie do chwili \(T\), liczba kroków jest w przybliżeniu równa:

\[ N=\frac{T-t_0}{h}. \]

Im mniejszy jest krok \(h\), tym gęściej rozmieszczone są punkty rozwiązania i zazwyczaj tym większa jest dokładność. Jednocześnie rośnie jednak liczba działań, koszt obliczeń oraz wpływ błędów zaokrągleń.

W metodach adaptacyjnych krok nie jest stały. Zamiast \(h\) stosuje się kolejne wartości:

\[ h_0,\ h_1,\ h_2,\ldots, \]

dostosowywane do lokalnego zachowania rozwiązania.

Całkowa interpretacja jednego kroku

Z równania:

\[ y'(t)=f(t,y(t)) \]

wynika:

\[ y(t_{n+1})-y(t_n) = \int_{t_n}^{t_{n+1}} f(t,y(t))\,dt. \]

A zatem:

\[ y(t_{n+1}) = y(t_n) + \int_{t_n}^{t_{n+1}} f(t,y(t))\,dt. \]

Wiele metod rozwiązywania równań różniczkowych można rozumieć jako różne sposoby przybliżania tej całki.

• Euler jawny wykorzystuje wartość funkcji na początku kroku.
• Euler niejawny wykorzystuje wartość na końcu kroku.
• Metoda trapezowa uśrednia wartości na początku i końcu.
• Metody Rungego – Kutty obliczają kilka pośrednich ocen nachylenia.

Związek ten łączy numeryczne rozwiązywanie równań różniczkowych z całkowaniem numerycznym.

Pole kierunków

W każdym punkcie płaszczyzny \((t,y)\) wartość:

\[ f(t,y) \]

określa nachylenie krzywej rozwiązania. Jeżeli w wielu punktach narysujemy krótkie odcinki o takim nachyleniu, otrzymamy pole kierunków.

Krzywe rozwiązań powinny w każdym swoim punkcie przebiegać zgodnie z lokalnym kierunkiem pola. Metody numeryczne można interpretować jako wędrówkę po tym polu: od aktualnego punktu przechodzimy do następnego, korzystając z jednej lub kilku ocen kierunku.

Jawna metoda Eulera

Jawna metoda Eulera jest najprostszą metodą jednokrokową. W aktualnym punkcie:

\[ (t_n,y_n) \]

obliczamy nachylenie:

\[ f(t_n,y_n). \]

Następnie przesuwamy się o krok \(h\) wzdłuż prostej o tym nachyleniu:

\[ y_{n+1} = y_n+h f(t_n,y_n). \]

Metodę tę można otrzymać z rozwinięcia Taylora:

\[ y(t_n+h) = y(t_n) + hy'(t_n) + \frac{h^2}{2}y''(\xi_n), \]

gdzie \(\xi_n\in(t_n,t_{n+1})\). Po pominięciu składnika zawierającego \(h^2\) oraz wykorzystaniu równania \(y'=f(t,y)\) otrzymujemy wzór Eulera.

Algorytm jawnej metody Eulera

1. Ustal punkt początkowy \((t_0,y_0)\).
2. Wybierz długość kroku \(h\).
3. Oblicz nachylenie \(f(t_n,y_n)\).
4. Wyznacz:
   \[
   y_{n+1}=y_n+h f(t_n,y_n).
   \]
5. Przyjmij:
   \[
   t_{n+1}=t_n+h.
   \]
6. Powtarzaj obliczenia do osiągnięcia końca przedziału.

Przykład jawnej metody Eulera

Rozważmy zagadnienie:

\[ y'=-2y, \qquad y(0)=1. \]

Przyjmujemy krok:

\[ h=0{,}25. \]

Wzór Eulera przyjmuje postać:

\[ y_{n+1} = y_n+0{,}25(-2y_n) = 0{,}5y_n. \]

Kolejne wartości wynoszą więc:

\[ y_0=1, \qquad y_1=0{,}5, \qquad y_2=0{,}25, \qquad y_3=0{,}125, \qquad y_4=0{,}0625. \]

Dokładna wartość w chwili \(t=1\) wynosi:

\[ y(1)=e^{-2}\approx0{,}1353352832. \]

Euler jawny daje:

\[ y_4=0{,}0625, \]

więc błąd jest wyraźnie widoczny. Zmniejszenie kroku poprawiłoby dokładność.

Dokładność jawnej metody Eulera

Lokalny błąd obcięcia, czyli błąd powstający podczas jednego kroku rozpoczętego od dokładnej wartości, jest rzędu:

\[ O(h^2). \]

Błędy poszczególnych kroków kumulują się. Globalny błąd rozwiązania na skończonym przedziale jest rzędu:

\[ O(h). \]

Dwukrotne zmniejszenie kroku powinno więc w typowej sytuacji zmniejszyć główny składnik błędu globalnego około dwa razy.

Niejawna metoda Eulera

W niejawnej metodzie Eulera wykorzystujemy nachylenie w punkcie końcowym kroku:

\[ y_{n+1} = y_n+ h f(t_{n+1},y_{n+1}). \]

Nowa wartość \(y_{n+1}\) występuje po obu stronach równania. Metoda nie podaje więc wyniku bezpośrednio i w każdym kroku trzeba rozwiązać równanie algebraiczne.

Dla funkcji liniowej może być to proste. Dla równania:

\[ y'=-2y \]

otrzymujemy:

\[ y_{n+1} = y_n-2h y_{n+1}. \]

Stąd:

\[ (1+2h)y_{n+1}=y_n, \]

czyli:

\[ y_{n+1} = \frac{y_n}{1+2h}. \]

Dla \(h=0{,}25\):

\[ y_{n+1} = \frac23y_n. \]

Kolejne przybliżenia wynoszą:

\[ 1,\quad 0{,}6666666667,\quad 0{,}4444444444,\quad 0{,}2962962963,\quad 0{,}1975308642. \]

Niejawna metoda Eulera ma, podobnie jak wariant jawny, globalny błąd rzędu \(O(h)\). Jej najważniejszą zaletą nie jest wyższy rząd dokładności, lecz znacznie lepsza stabilność.

Rozwiązywanie równania wewnątrz kroku

Dla nieliniowej funkcji \(f\) zależność:

\[ y_{n+1} - y_n - h f(t_{n+1},y_{n+1}) = 0 \]

jest równaniem nieliniowym względem \(y_{n+1}\). Można je rozwiązywać metodą Newtona lub innym algorytmem opisanym w artykule Numeryczne rozwiązywanie równań nieliniowych.

W takim przypadku pojedynczy krok metody niejawnej zawiera wewnętrzną iterację z własnym kryterium stopu. Jest to istotna różnica w porównaniu z metodą jawną.

Ulepszona metoda Eulera – metoda Heuna

Jawna metoda Eulera korzysta tylko z nachylenia na początku kroku. Metoda Heuna oblicza dodatkowo przewidywane nachylenie na jego końcu, a następnie uśrednia oba kierunki.

Jest to jawny schemat typu predyktor – korektor.

Predykcja

Najpierw wykonujemy zwykły krok Eulera:

\[ \widetilde y_{n+1} = y_n+ h f(t_n,y_n). \]

Wartość \(\widetilde y_{n+1}\) jest przewidywanym położeniem na końcu kroku.

Korekcja

Obliczamy nachylenie w przewidzianym punkcie:

\[ f(t_{n+1},\widetilde y_{n+1}), \]

a następnie stosujemy średnią nachyleń:

\[ y_{n+1} = y_n+ \frac h2 \left[ f(t_n,y_n) + f(t_{n+1},\widetilde y_{n+1}) \right]. \]

Przykład metody Heuna

Dla równania:

\[ y'=-2y, \qquad y_0=1, \qquad h=0{,}25 \]

predykcja daje:

\[ \widetilde y_1 = 1+0{,}25(-2) = 0{,}5. \]

Nachylenia na początku i w punkcie przewidywanym wynoszą:

\[ k_1=-2, \qquad k_2=-2\cdot0{,}5=-1. \]

Po korekcji:

\[ y_1 = 1+ \frac{0{,}25}{2} (-2-1) = 0{,}625. \]

Kolejne przybliżenia są równe:

\[ 1,\quad 0{,}625,\quad 0{,}390625,\quad 0{,}244140625,\quad 0{,}1525878906. \]

Metoda Heuna ma lokalny błąd rzędu:

\[ O(h^3), \]

oraz globalny błąd rzędu:

\[ O(h^2). \]

Jawna metoda punktu środkowego

Inną metodą drugiego rzędu jest jawna metoda punktu środkowego. Najpierw obliczamy przybliżenie rozwiązania w połowie kroku:

\[ y_{n+\frac12} = y_n+ \frac h2f(t_n,y_n). \]

Następnie wykorzystujemy nachylenie w tym punkcie:

\[ y_{n+1} = y_n+ h f\left( t_n+\frac h2, y_{n+\frac12} \right). \]

Metoda Heuna i metoda punktu środkowego są różnymi metodami Rungego – Kutty drugiego rzędu. W obu wykonuje się dwie oceny funkcji w jednym kroku, lecz w innych punktach.

Metody Rungego – Kutty

Metody Rungego – Kutty zwiększają dokładność bez jawnego obliczania wyższych pochodnych rozwiązania. W jednym kroku wykonują kilka ocen funkcji \(f(t,y)\), a następnie łączą je w odpowiednio dobraną średnią ważoną.

Metoda Eulera może być traktowana jako metoda Rungego – Kutty pierwszego rzędu. Heun i metoda punktu środkowego należą do metod drugiego rzędu.

Klasyczna metoda Rungego – Kutty czwartego rzędu – RK4

Najbardziej znanym schematem jest klasyczna metoda Rungego – Kutty czwartego rzędu, oznaczana jako RK4.

W każdym kroku obliczamy cztery nachylenia.

Nachylenie na początku kroku:

\[ k_1=f(t_n,y_n). \]

Pierwsza ocena w połowie kroku:

\[ k_2= f\left( t_n+\frac h2, y_n+\frac h2k_1 \right). \]

Druga ocena w połowie kroku:

\[ k_3= f\left( t_n+\frac h2, y_n+\frac h2k_2 \right). \]

Nachylenie na końcu kroku:

\[ k_4= f(t_n+h,y_n+hk_3). \]

Nową wartość obliczamy ze wzoru:

\[ y_{n+1} = y_n+ \frac h6 \left( k_1+2k_2+2k_3+k_4 \right). \]

Nachylenia \(k_2\) i \(k_3\), obliczane w środku kroku, otrzymują dwukrotnie większe wagi niż wartości na początku i końcu.

Przykład jednego kroku RK4

Dla równania:

\[ y'=-2y, \qquad y_0=1, \qquad h=0{,}25 \]

otrzymujemy:

\[ k_1=-2, \]

\[ k_2= -2\left(1+\frac{0{,}25}{2}(-2)\right) =-1{,}5, \]

\[ k_3= -2\left(1+\frac{0{,}25}{2}(-1{,}5)\right) =-1{,}625, \]

\[ k_4= -2\left(1+0{,}25(-1{,}625)\right) =-1{,}1875. \]

Stąd:

\[ y_1 = 1+ \frac{0{,}25}{6} \left[ -2+2(-1{,}5)+2(-1{,}625)-1{,}1875 \right] \]

\[ y_1 \approx 0{,}6067708333. \]

Dokładna wartość wynosi:

\[ y(0{,}25)=e^{-0{,}5} \approx 0{,}6065306597. \]

Już po jednym kroku wynik jest bardzo dokładny.

Dokładność RK4

Lokalny błąd obcięcia klasycznej metody RK4 jest rzędu:

\[ O(h^5), \]

a błąd globalny:

\[ O(h^4). \]

Dwukrotne zmniejszenie kroku powinno więc w typowej sytuacji zmniejszyć główny składnik błędu globalnego około:

\[ 2^4=16 \]

razy.

Porównanie metod na wspólnym przykładzie

Porównajmy wyniki dla zagadnienia:

\[ y'=-2y, \qquad y(0)=1, \]

przy kroku:

\[ h=0{,}25. \]

Błędy bezwzględne w chwili \(t=1\) wynoszą w przybliżeniu:

RK4 wykonuje więcej obliczeń w jednym kroku, ale uzyskuje znacznie większą dokładność.

Błąd lokalny i globalny

Przy analizie metod rozwiązywania równań różniczkowych należy rozróżnić dwa podstawowe rodzaje błędów.

Lokalny błąd obcięcia

Jest to błąd wykonany w pojedynczym kroku przy założeniu, że rozpoczynamy go od dokładnej wartości:

\[ y_n=y(t_n). \]

Błąd globalny

Jest to różnica pomiędzy dokładnym rozwiązaniem i wartością otrzymaną po wielu krokach:

\[ e_n=y(t_n)-y_n. \]

Błąd globalny zawiera wpływ błędów powstałych we wszystkich wcześniejszych krokach.

Czy istnieje kryterium stopu?

W metodach rozwiązywania równania:

\[ f(x)=0 \]

kolejne iteracje prowadzi się zazwyczaj do chwili, gdy przybliżenie pierwiastka spełni zadane kryterium dokładności.

Przy numerycznym rozwiązywaniu równania różniczkowego sytuacja jest inna. Najczęściej z góry ustalamy przedział:

\[ [t_0,T] \]

i obliczamy rozwiązanie aż do osiągnięcia chwili końcowej:

\[ t_N=T. \]

Nie próbujemy uzyskać jednego punktu granicznego przez nieskończone poprawianie przybliżenia. Każdy krok wyznacza kolejny fragment trajektorii.

Obliczenia mogą zakończyć się również wtedy, gdy:

• rozwiązanie osiągnie określoną wartość,
• zajdzie zdarzenie, na przykład \(y(t)=0\),
• trajektoria opuści dopuszczalny obszar,
• zostanie osiągnięta maksymalna liczba kroków,
• rozwiązanie zbliży się do stanu ustalonego,
• długość kroku spadnie poniżej dopuszczalnej wartości.

W metodach adaptacyjnych tolerancja błędu służy przede wszystkim do przyjmowania lub odrzucania pojedynczego kroku oraz doboru jego długości, a nie do określania końca całej trajektorii.

Metody adaptacyjne Rungego – Kutty

Stały krok nie zawsze jest efektywny. Na fragmentach, gdzie rozwiązanie zmienia się powoli, można stosować duży krok. W miejscach gwałtownej zmiany potrzebny jest krok mniejszy.

Popularne metody adaptacyjne wykorzystują dwie formuły Rungego – Kutty różnych rzędów, obliczane przy użyciu wspólnych wartości pośrednich. Przykładami są:

• metoda Rungego – Kutty – Fehlberga,
• metoda Dormanda – Prince’a,
• schematy RK45.

Jeżeli w jednym kroku otrzymano dwa przybliżenia:

\[ y_{n+1}^{(p)} \qquad\text{oraz}\qquad y_{n+1}^{(p+1)}, \]

ich różnica służy jako oszacowanie błędu:

\[ e_n \approx \left| y_{n+1}^{(p+1)} - y_{n+1}^{(p)} \right|. \]

Jeżeli błąd przekracza tolerancję, krok zostaje odrzucony i wykonany ponownie z mniejszym \(h\). Jeżeli jest znacznie mniejszy od tolerancji, kolejny krok można zwiększyć.

W efekcie punkty rozwiązania mogą być rozmieszczone nierównomiernie:

\[ t_{n+1}=t_n+h_n. \]

Stabilność metody numerycznej

Wysoki rząd dokładności nie gwarantuje poprawnego zachowania metody dla dowolnego kroku. Algorytm może być dokładny dla małego \(h\), ale niestabilny po jego zwiększeniu.

Do analizy stabilności stosuje się równanie testowe:

\[ y'=\lambda y, \qquad \operatorname{Re}(\lambda)<0. \]

Rozwiązanie dokładne:

\[ y(t)=y_0e^{\lambda t} \]

zanika do zera.

Stabilność jawnej metody Eulera

Euler jawny daje:

\[ y_{n+1} = (1+h\lambda)y_n. \]

Aby kolejne wartości zanikały, potrzebny jest warunek:

\[ |1+h\lambda|<1. \]

Dla równania:

\[ y'=-15y \]

otrzymujemy:

\[ |1-15h|<1, \]

czyli:

\[ 0<h<\frac{2}{15} \approx0{,}1333. \]

Dla \(h=0{,}2\):

\[ 1-15h=-2. \]

Kolejne wartości wynoszą:

\[ 1,\ -2,\ 4,\ -8,\ 16,\ldots \]

Metoda oscyluje i rozbiega się, mimo że rozwiązanie dokładne bardzo szybko zanika.

Stabilność niejawnej metody Eulera

Euler niejawny daje:

\[ y_{n+1} = \frac{y_n}{1-h\lambda}. \]

Dla \(\lambda=-15\) oraz \(h=0{,}2\):

\[ y_{n+1} = \frac{y_n}{4}. \]

Kolejne wartości wynoszą:

\[ 1,\quad 0{,}25,\quad 0{,}0625,\quad 0{,}015625,\ldots \]

Rozwiązanie numeryczne zachowuje stabilny zanik.

Stabilność nie jest tym samym co dokładność. Metoda niejawna może zachowywać właściwy charakter rozwiązania przy dużym kroku, ale jej wartości nadal mogą wyraźnie różnić się od rozwiązania dokładnego.

Równania sztywne

Równanie lub układ równań nazywa się sztywnym, gdy w rozwiązaniu występują procesy zachodzące w bardzo różnych skalach czasowych. Niektóre składniki zanikają bardzo szybko, podczas gdy inne zmieniają się powoli.

Jawna metoda może wówczas wymagać niezwykle małego kroku wyłącznie ze względu na stabilność, mimo że dokładny przebieg rozwiązania jest już prawie gładki.

W takich problemach często stosuje się:

• niejawną metodę Eulera,
• niejawną metodę trapezową,
• metody Adamsa – Moultona,
• wsteczne formuły różnicowe BDF,
• niejawne metody Rungego – Kutty.

Niejawna metoda trapezowa

Metoda trapezowa wykorzystuje średnią nachyleń na początku i końcu kroku:

\[ y_{n+1} = y_n+ \frac h2 \left[ f(t_n,y_n) + f(t_{n+1},y_{n+1}) \right]. \]

W przeciwieństwie do metody Heuna wartość \(y_{n+1}\) wewnątrz funkcji nie jest tylko predykcją, lecz poszukiwanym wynikiem. Jest to więc metoda niejawna, wymagająca rozwiązania równania w każdym kroku.

Metoda trapezowa ma globalny błąd rzędu \(O(h^2)\) i bardzo dobre właściwości stabilnościowe.

Metody wielokrokowe

Metody Eulera i Rungego – Kutty są metodami jednokrokowymi: do wyznaczenia \(y_{n+1}\) wykorzystują przede wszystkim informacje z bieżącego kroku.

Metody wielokrokowe korzystają również z wcześniej obliczonych wartości. Pozwala to ograniczyć liczbę nowych obliczeń funkcji \(f\).

Jawna metoda Adamsa – Bashfortha

Dwukrokowa metoda Adamsa – Bashfortha ma postać:

\[ y_{n+1} = y_n+ \frac h2 \left[ 3f(t_n,y_n) - f(t_{n-1},y_{n-1}) \right]. \]

Do rozpoczęcia obliczeń potrzebna jest dodatkowa wartość \(y_1\), wyznaczona na przykład metodą Rungego – Kutty.

Metody Adamsa – Moultona

Metody Adamsa – Moultona są schematami niejawnymi. Wartość funkcji w punkcie końcowym zależy w nich od poszukiwanego \(y_{n+1}\).

Często łączy się metodę Adamsa – Bashfortha jako predyktor z metodą Adamsa – Moultona jako korektor.

Wsteczne formuły różnicowe BDF

Metody BDF przybliżają pochodną za pomocą wartości rozwiązania z kilku kolejnych chwil. Dwukrokowy wzór BDF ma postać:

\[ \frac{ 3y_{n+1}-4y_n+y_{n-1} }{ 2h } = f(t_{n+1},y_{n+1}). \]

Metody te są niejawne i szczególnie użyteczne w rozwiązywaniu układów sztywnych.

Układy równań różniczkowych

Większość metod można bezpośrednio stosować do układów równań. Zamiast jednej liczby \(y\) mamy wektor:

\[ \mathbf y(t) = \begin{bmatrix} y_1(t)\\ y_2(t)\\ \vdots\\ y_m(t) \end{bmatrix}. \]

Układ zapisujemy jako:

\[ \mathbf y'(t) = \mathbf f(t,\mathbf y). \]

Jawna metoda Eulera przyjmuje wtedy postać:

\[ \mathbf y_{n+1} = \mathbf y_n+ h\mathbf f(t_n,\mathbf y_n). \]

Wzory Rungego – Kutty pozostają analogiczne, lecz wartości \(k_i\) są wektorami.

Równania wyższego rzędu

Równanie drugiego rzędu:

\[ y''=g(t,y,y') \]

można sprowadzić do układu dwóch równań pierwszego rzędu. Wprowadzamy:

\[ y_1=y, \qquad y_2=y'. \]

Wtedy:

\[ y_1'=y_2, \]

\[ y_2'=g(t,y_1,y_2). \]

Dzięki temu tymi samymi metodami można rozwiązywać między innymi równania ruchu mechanicznego:

\[ m y''+c y'+ky=F(t). \]

Zagadnienia brzegowe

Artykuł koncentruje się na zagadnieniach początkowych, w których wszystkie warunki są podane w jednym punkcie \(t_0\).

W zagadnieniach brzegowych warunki są zadane w różnych punktach, na przykład:

\[ y(a)=\alpha, \qquad y(b)=\beta. \]

Do ich rozwiązywania stosuje się między innymi:

• metodę strzałów,
• metody różnic skończonych,
• metody elementów skończonych,
• metody kolokacji.

W metodzie strzałów zagadnienie brzegowe sprowadza się do wielokrotnego rozwiązywania zagadnienia początkowego i dobierania brakującego warunku tak, aby spełnić warunek na drugim końcu przedziału.

Wpływ błędów zaokrągleń

Zmniejszanie kroku redukuje błąd obcięcia, ale zwiększa liczbę działań. Każde działanie jest wykonywane w skończonej precyzji, dlatego bardzo duża liczba kroków może powodować kumulowanie błędów zaokrągleń.

Znaczenie ograniczonej reprezentacji liczb zostało omówione w artykule Arytmetyka zmiennoprzecinkowa – cecha, mantysa i zapis binarny.

W praktyce krok powinien być wystarczająco mały, aby ograniczyć błąd metody, ale nie mniejszy, niż wymaga tego oczekiwana dokładność.

Porównanie najważniejszych metod

Którą metodę wybrać?

• Euler jawny jest odpowiedni do nauki, prostych symulacji i wstępnego sprawdzania modeli.
• Euler niejawny jest użyteczny, gdy najważniejsza jest stabilność, szczególnie dla równań sztywnych.
• Heun i metoda punktu środkowego zapewniają dobrą równowagę między prostotą i dokładnością.
• RK4 jest dobrym wyborem dla wielu niesztwnych zagadnień z ustalonym krokiem.
• RK45 i inne metody adaptacyjne sprawdzają się, gdy rozwiązanie zmienia się z różną szybkością.
• Metody BDF są często stosowane w dużych, sztywnych układach.

Praktyczne zasady obliczeń

• Sprawdź, czy rozwiązujesz zagadnienie początkowe, czy brzegowe.
• Ustal przedział obliczeń i oczekiwaną dokładność.
• Dobierz krok odpowiedni do szybkości zmian rozwiązania.
• Porównuj wyniki dla kroków \(h\) i \(h/2\).
• Nie utożsamiaj stabilności z dokładnością.
• Dla metod niejawnych kontroluj dokładność wewnętrznego rozwiązania równania.
• W metodach adaptacyjnych ustal tolerancję bezwzględną i względną.
• Kontroluj, czy rozwiązanie nie opuszcza dziedziny modelu.
• Uwzględniaj możliwość wystąpienia zdarzeń pomiędzy punktami siatki.
• Przy układach sztywnych nie zwiększaj bezkrytycznie kroku metody jawnej.
• Sprawdzaj zachowanie wielkości, które powinny być zachowane, na przykład energii lub masy.
• Porównuj wynik z rozwiązaniem dokładnym lub inną metodą, jeżeli jest to możliwe.

Numeryczne rozwiązanie równania różniczkowego nie jest zazwyczaj wzorem funkcji, lecz obliczoną krok po kroku trajektorią jej przybliżonych wartości.

Podsumowanie

Numeryczne rozwiązywanie równań różniczkowych polega najczęściej na wyznaczaniu kolejnych wartości rozwiązania w punktach siatki czasowej. Zamiast jawnej postaci funkcji otrzymujemy dyskretny zbiór punktów opisujących jej przybliżoną trajektorię.

Jawna metoda Eulera wykorzystuje nachylenie na początku kroku. Jest bardzo prosta, lecz ma dokładność pierwszego rzędu i może być niestabilna.

Niejawna metoda Eulera korzysta z nachylenia w punkcie końcowym. Wymaga rozwiązania równania w każdym kroku, ale ma znacznie lepsze właściwości stabilnościowe.

Metoda Heuna i metoda punktu środkowego wykonują po dwie oceny nachylenia i osiągają dokładność drugiego rzędu.

Klasyczna metoda RK4 wykorzystuje cztery nachylenia w jednym kroku. Jej globalny błąd jest rzędu \(h^4\), dzięki czemu zapewnia dużą dokładność bez konieczności obliczania wyższych pochodnych.

Metody adaptacyjne samodzielnie zmieniają długość kroku na podstawie oszacowania lokalnego błędu. Metody niejawne i wielokrokowe są szczególnie istotne w rozwiązywaniu równań sztywnych.

Obliczenia prowadzi się zwykle do końca ustalonego przedziału albo do wystąpienia określonego zdarzenia. Klasyczne kryterium stopu pojawia się głównie wewnątrz metod niejawnych, gdy w pojedynczym kroku trzeba iteracyjnie rozwiązać równanie algebraiczne.

Utworzono: 17.06.2026 | Zmodyfikowano: 26.07.2026

Całkowanie numeryczne – metoda prostokątów, trapezów i Simpsona

Mgr inż.

Sebastian Dziarmaga-Działyński

Wirtualna kancelaria korepetytorska i konsultacyjna – wszechwiedza.pl

Całka oznaczona może być interpretowana jako zorientowane pole obszaru pomiędzy wykresem funkcji a osią argumentów. Jeżeli znamy funkcję pierwotną, jej wartość obliczamy za pomocą twierdzenia Newtona – Leibniza:

\[ \int_a^b f(x)\,dx=F(b)-F(a), \qquad F'(x)=f(x). \]

Nie każda funkcja ma jednak funkcję pierwotną, którą można wyrazić za pomocą funkcji elementarnych. Ponadto w praktycznych zastosowaniach funkcja może być znana jedynie w postaci tabeli wyników pomiarów albo jej obliczanie może być bardzo kosztowne.

W takich sytuacjach stosuje się całkowanie numeryczne, nazywane również kwadraturą numeryczną. Pole pod wykresem zastępuje się wówczas sumą pól prostszych figur: prostokątów, trapezów lub obszarów ograniczonych parabolami interpolacyjnymi.

Całkowanie numeryczne jest jednym z podstawowych zagadnień metod numerycznych. Jego dokładność zależy od zastosowanej kwadratury, liczby podprzedziałów, własności funkcji oraz błędów obliczeń komputerowych.

Całka bez elementarnej funkcji pierwotnej

Prostym przykładem pokazującym użyteczność metod numerycznych jest całka:

\[ \int_0^1 e^{x^2}\,dx. \]

Funkcja \(e^{x^2}\) nie ma funkcji pierwotnej wyrażalnej za pomocą funkcji elementarnych. Można ją zapisać za pomocą funkcji specjalnej \(\operatorname{erfi}(x)\), ale do wyznaczenia liczbowej wartości całki oznaczonej nie jest to konieczne.

Metody numeryczne pozwalają otrzymać:

\[ \int_0^1 e^{x^2}\,dx \approx 1{,}4626517459. \]

Nie musimy więc znać elementarnego wzoru na funkcję pierwotną, aby z dużą dokładnością obliczyć wartość całki oznaczonej.

Podział przedziału całkowania

Rozważmy całkę:

\[ I=\int_a^b f(x)\,dx. \]

Przedział \([a,b]\) dzielimy na \(n\) równych części. Długość każdego podprzedziału wynosi:

\[ h=\frac{b-a}{n}. \]

Węzły podziału mają współrzędne:

\[ x_i=a+ih, \qquad i=0,1,\ldots,n. \]

Wartości funkcji w węzłach będziemy oznaczać jako:

\[ f_i=f(x_i). \]

Metody prostokątów, trapezów i Simpsona różnią się sposobem zastępowania funkcji na każdym fragmencie przedziału.

Błąd kwadratury numerycznej

Jeżeli \(Q_n\) oznacza wynik otrzymany za pomocą kwadratury numerycznej, jej błąd można zapisać jako:

\[ E_n=I-Q_n. \]

Najczęściej interesuje nas moduł błędu:

\[ |E_n|=|I-Q_n|. \]

Wzory na błędy metod całkowania zawierają największą wartość modułu odpowiedniej pochodnej na całym przedziale. Wprowadzamy oznaczenie:

\[ M_k= \sup_{x\in[a,b]} \left|f^{(k)}(x)\right|. \]

Supremum jest najmniejszym ograniczeniem górnym zbioru wartości. Jeżeli pochodna \(f^{(k)}\) jest ciągła na domkniętym przedziale \([a,b]\), supremum jest osiągane i można po prostu napisać:

\[ M_k= \max_{x\in[a,b]} \left|f^{(k)}(x)\right|. \]

Nie sprawdzamy więc pochodnej tylko w kilku węzłach obliczeniowych. Poszukujemy największego modułu pochodnej wśród wszystkich punktów przedziału całkowania.

Oszacowanie błędu jest ograniczeniem z góry. Gwarantuje, że błąd nie przekroczy określonej wartości, ale rzeczywisty błąd może być — i zazwyczaj jest — znacznie mniejszy.

Metoda prostokątów

W metodzie prostokątów funkcję na każdym podprzedziale zastępuje się funkcją stałą. Pole pod wykresem przybliża się sumą pól prostokątów o szerokości \(h\).

W zależności od miejsca, w którym obliczamy wysokość prostokąta, wyróżniamy metodę prostokątów lewych, prawych i środkowych.

Prostokąty lewe

Wysokość prostokąta na podprzedziale \([x_i,x_{i+1}]\) jest równa wartości funkcji na jego lewym końcu:

\[ f(x_i). \]

Złożony wzór lewych prostokątów ma postać:

\[ L_n= h\sum_{i=0}^{n-1}f(x_i). \]

Jeżeli funkcja jest malejąca, lewe prostokąty zwykle zawyżają pole. Dla funkcji rosnącej zazwyczaj je zaniżają.

Prostokąty prawe

Wysokość prostokąta wyznacza wartość funkcji na prawym końcu podprzedziału:

\[ f(x_{i+1}). \]

Złożony wzór prawych prostokątów jest równy:

\[ P_n= h\sum_{i=1}^{n}f(x_i). \]

Dla funkcji malejącej prawe prostokąty zazwyczaj zaniżają pole, natomiast dla funkcji rosnącej je zawyżają.

Prostokąty środkowe

W metodzie punktu środkowego wysokość prostokąta obliczamy w środku każdego podprzedziału:

\[ m_i=\frac{x_i+x_{i+1}}{2} =x_i+\frac{h}{2}. \]

Otrzymujemy wzór:

\[ M_n= h\sum_{i=0}^{n-1} f\left(x_i+\frac{h}{2}\right). \]

Metoda środkowych prostokątów jest zazwyczaj znacznie dokładniejsza od wariantu lewego i prawego. Wynika to między innymi z częściowego znoszenia się błędów po obu stronach środka podprzedziału.

Oszacowanie błędu metody prostokątów

Dla lewych i prawych prostokątów:

\[ |E_{L,P}| \leq \frac{b-a}{2}hM_1. \]

Po podstawieniu \(h=(b-a)/n\):

\[ |E_{L,P}| \leq \frac{(b-a)^2}{2n}M_1. \]

Błąd maleje proporcjonalnie do \(h\), dlatego są to metody rzędu:

\[ O(h). \]

Dla prostokątów środkowych:

\[ |E_M| \leq \frac{b-a}{24}h^2M_2, \]

czyli:

\[ |E_M| \leq \frac{(b-a)^3}{24n^2}M_2. \]

Metoda punktu środkowego ma dokładność rzędu:

\[ O(h^2). \]

Metoda trapezów

W metodzie trapezów funkcję na każdym podprzedziale zastępuje się odcinkiem łączącym punkty:

\[ (x_i,f(x_i)) \qquad\text{oraz}\qquad (x_{i+1},f(x_{i+1})). \]

Jest to interpolacja wielomianem pierwszego stopnia. Obszar pod odcinkiem ma kształt trapezu.

Pole pojedynczego trapezu wynosi:

\[ T_i= h\frac{f(x_i)+f(x_{i+1})}{2}. \]

Po zsumowaniu pól wszystkich trapezów otrzymujemy złożony wzór:

\[ T_n= \frac{h}{2} \left[ f(x_0) + 2\sum_{i=1}^{n-1}f(x_i) + f(x_n) \right]. \]

Trapez a prostokąt środkowy

Metoda trapezów może na pierwszy rzut oka przypominać metodę prostokątów środkowych. W obu przypadkach można bowiem mówić o pewnej „średniej wysokości” figury. Są to jednak dwie różne wielkości.

W metodzie prostokątów środkowych wysokość jest wartością funkcji dla argumentu znajdującego się w środku podprzedziału:

\[ H_M= f\left(\frac{x_i+x_{i+1}}{2}\right). \]

W metodzie trapezów średnia wysokość jest średnią arytmetyczną wartości funkcji na obu końcach podprzedziału:

\[ H_T= \frac{f(x_i)+f(x_{i+1})}{2}. \]

W ogólnym przypadku:

\[ f\left(\frac{x_i+x_{i+1}}{2}\right) \ne \frac{f(x_i)+f(x_{i+1})}{2}. \]

Równość zachodzi dla funkcji liniowych, ale nie musi zachodzić dla funkcji nieliniowej.

W prostokącie środkowym najpierw uśredniamy argumenty, a następnie obliczamy wartość funkcji. W trapezie najpierw obliczamy wartości funkcji na końcach, a następnie je uśredniamy.

Jeżeli funkcja jest wypukła na podprzedziale, zachodzi nierówność:

\[ f\left(\frac{x_i+x_{i+1}}{2}\right) \leq \frac{f(x_i)+f(x_{i+1})}{2}. \]

Wtedy trapez ma większą średnią wysokość niż prostokąt środkowy. Dla funkcji wklęsłej kierunek nierówności jest przeciwny.

Błąd metody trapezów

Jeżeli funkcja ma ciągłą drugą pochodną, błąd złożonej metody trapezów spełnia oszacowanie:

\[ |E_T| \leq \frac{b-a}{12}h^2M_2. \]

Równoważnie:

\[ |E_T| \leq \frac{(b-a)^3}{12n^2}M_2. \]

Metoda trapezów, podobnie jak metoda prostokątów środkowych, ma błąd rzędu:

\[ O(h^2). \]

Współczynnik w górnym oszacowaniu błędu metody punktu środkowego jest jednak dwukrotnie mniejszy niż dla trapezów:

\[ \frac{1}{24} \frac{1}{12}. \]

Nie oznacza to, że w każdym pojedynczym przykładzie prostokąty środkowe muszą dać dokładnie dwukrotnie mniejszy błąd. Są to oszacowania najgorszego możliwego przypadku.

Metoda Simpsona

Metoda Simpsona zastępuje funkcję na każdej parze sąsiednich podprzedziałów wielomianem drugiego stopnia. Parabola interpolacyjna przechodzi dokładnie przez trzy kolejne węzły:

\[ (x_{2k},f_{2k}), \qquad (x_{2k+1},f_{2k+1}), \qquad (x_{2k+2},f_{2k+2}). \]

Interpolanta nie musi być styczna do wykresu całkowanej funkcji. Ma jedynie przechodzić przez trzy wybrane punkty. Pomiędzy węzłami może znajdować się zarówno nad wykresem, jak i pod nim oraz przecinać go w węźle środkowym.

Dla trzech kolejnych węzłów \(x_0\), \(x_1\), \(x_2\), oddalonych o \(h\), pojedynczy wzór Simpsona ma postać:

\[ \int_{x_0}^{x_2}f(x)\,dx \approx \frac{h}{3} \left[ f(x_0)+4f(x_1)+f(x_2) \right]. \]

Aby można było grupować podprzedziały parami, liczba \(n\) musi być parzysta.

Złożony wzór Simpsona wynosi:

\[ S_n= \frac{h}{3} \left[ f(x_0) + 4\sum_{\substack{i=1\\ i\ \mathrm{nieparzyste}}}^{n-1}f(x_i) + 2\sum_{\substack{i=2\\ i\ \mathrm{parzyste}}}^{n-2}f(x_i) + f(x_n) \right]. \]

Wagi kolejnych wartości funkcji mają więc układ:

\[ 1,\ 4,\ 2,\ 4,\ 2,\ldots,\ 4,\ 1. \]

Dlaczego Simpson jest tak dokładny?

Choć konstrukcja metody opiera się na interpolacji wielomianem drugiego stopnia, wzór Simpsona jest dokładny również dla wielomianów stopnia trzeciego. Wynika to z symetrii węzłów i całkowania na parze podprzedziałów.

Jeżeli funkcja ma ciągłą czwartą pochodną, błąd złożonej metody Simpsona spełnia:

\[ |E_S| \leq \frac{b-a}{180}h^4M_4. \]

Po zapisaniu przez liczbę podprzedziałów:

\[ |E_S| \leq \frac{(b-a)^5}{180n^4}M_4. \]

Metoda Simpsona ma zatem dokładność rzędu:

\[ O(h^4). \]

Dwukrotne zmniejszenie kroku powoduje w typowej sytuacji zmniejszenie głównego składnika błędu około:

\[ 2^4=16 \]

razy.

Porównanie metod na wspólnym przykładzie

Rozważmy całkę:

\[ I= \int_0^1\frac{1}{1+x^2}\,dx. \]

W tym przypadku znamy wartość dokładną:

\[ I= \arctan 1-\arctan 0 = \frac{\pi}{4} \approx 0{,}7853981634. \]

Przyjmujemy cztery podprzedziały:

\[ n=4, \qquad h=\frac14. \]

Węzły są równe:

\[ x_0=0,\quad x_1=0{,}25,\quad x_2=0{,}5,\quad x_3=0{,}75,\quad x_4=1. \]

Dla funkcji:

\[ f(x)=\frac{1}{1+x^2} \]

otrzymujemy:

Już przy zaledwie czterech podprzedziałach metoda Simpsona daje wynik różniący się od wartości dokładnej o około sześć milionowych.

Czy oszacowania błędów są spełnione?

Dla funkcji:

\[ f(x)=\frac{1}{1+x^2} \]

potrzebne pochodne wynoszą:

\[ f'(x)= -\frac{2x}{(1+x^2)^2}, \]

\[ f''(x)= \frac{2(3x^2-1)}{(1+x^2)^3}, \]

\[ f^{(4)}(x)= \frac{24(5x^4-10x^2+1)} {(1+x^2)^5}. \]

Na przedziale \([0,1]\):

\[ M_1= \max_{x\in[0,1]}|f'(x)| = \frac{3\sqrt3}{8} \approx 0{,}649519, \]

\[ M_2= \max_{x\in[0,1]}|f''(x)| = 2, \]

\[ M_4= \max_{x\in[0,1]}|f^{(4)}(x)| = 24. \]

Dla \(h=1/4\) otrzymujemy następujące górne oszacowania:

W każdym przypadku rzeczywisty błąd jest mniejszy od wartości wynikającej z oszacowania. Widać też, że granice są dość zachowawcze, szczególnie dla metody Simpsona.

Nie jest to wada wzoru. Jego zadaniem jest zagwarantowanie maksymalnego błędu na podstawie najgorszej wartości odpowiedniej pochodnej na całym przedziale.

Dobór liczby podprzedziałów

Jeżeli chcemy zagwarantować błąd nie większy niż \(\varepsilon\), ze wzorów na błędy można wyznaczyć minimalną liczbę podprzedziałów.

Dla lewych i prawych prostokątów:

\[ n \geq \frac{(b-a)^2M_1}{2\varepsilon}. \]

Dla prostokątów środkowych:

\[ n \geq \sqrt{ \frac{(b-a)^3M_2}{24\varepsilon} }. \]

Dla trapezów:

\[ n \geq \sqrt{ \frac{(b-a)^3M_2}{12\varepsilon} }. \]

Dla metody Simpsona:

\[ n \geq \sqrt[4]{ \frac{(b-a)^5M_4}{180\varepsilon} }. \]

W metodzie Simpsona wynik zaokrąglamy w górę do najbliższej liczby parzystej.

Metoda Simpsona 3/8

Metoda Simpsona 3/8 jest kolejną zamkniętą kwadraturą Newtona – Cotesa. Wykorzystuje cztery kolejne węzły i wielomian interpolacyjny stopnia trzeciego.

Dla trzech podprzedziałów:

\[ \int_{x_0}^{x_3}f(x)\,dx \approx \frac{3h}{8} \left[ f(x_0)+3f(x_1)+3f(x_2)+f(x_3) \right]. \]

W wersji złożonej liczba podprzedziałów musi być podzielna przez trzy:

\[ Q_{3/8} = \frac{3h}{8} \left[ f(x_0)+f(x_n) + 3\sum_{\substack{1\leq i\leq n-1\\3\nmid i}}f(x_i) + 2\sum_{\substack{3\leq i\leq n-3\\3\mid i}}f(x_i) \right]. \]

Błąd można oszacować jako:

\[ |E_{3/8}| \leq \frac{b-a}{80}h^4M_4. \]

Metoda 3/8, podobnie jak klasyczny wzór Simpsona 1/3, ma dokładność rzędu \(O(h^4)\). Użycie interpolacji sześciennej nie powoduje automatycznie uzyskania wyższego rzędu dokładności. Klasyczna metoda Simpsona 1/3 ma nawet korzystniejszą stałą w oszacowaniu błędu.

Metoda Romberga

Metoda Romberga rozpoczyna obliczenia od kolejnych przybliżeń uzyskanych metodą trapezów dla kroków:

\[ h,\quad \frac{h}{2},\quad \frac{h}{4},\quad \frac{h}{8},\ldots \]

Następnie wykorzystuje ekstrapolację Richardsona do systematycznego usuwania dominujących składników błędu.

Nie chodzi więc jedynie o lepsze zaokrąglenie wyniku. Metoda wykorzystuje znaną strukturę błędu kwadratury trapezów i eliminuje kolejno składniki proporcjonalne do \(h^2\), \(h^4\), \(h^6\) i dalszych parzystych potęg kroku.

Pierwsza kolumna tablicy Romberga zawiera wyniki trapezów:

\[ R_{k,0}=T_{2^k}. \]

Kolejne elementy oblicza się ze wzoru:

\[ R_{k,j} = R_{k,j-1} + \frac{ R_{k,j-1}-R_{k-1,j-1} }{ 4^j-1 }. \]

Dla naszego przykładu:

\[ T_4=0{,}7827941176, \qquad T_8=0{,}7847471236. \]

Pierwsza ekstrapolacja daje:

\[ R= T_8+\frac{T_8-T_4}{3} = 0{,}7853981256. \]

Wartość dokładna wynosi \(0{,}7853981634\), więc już jedno poprawienie dwóch wyników metody trapezów daje bardzo dużą dokładność.

Kwadratury Gaussa

W klasycznych wzorach Newtona – Cotesa węzły są zwykle rozmieszczone równomiernie. W kwadraturach Gaussa zarówno położenie węzłów, jak i ich wagi dobiera się w sposób optymalny.

Dla kwadratury Gaussa – Legendre’a całkę sprowadza się najpierw do przedziału \([-1,1]\):

\[ x= \frac{a+b}{2} + \frac{b-a}{2}t. \]

Wtedy:

\[ \int_a^b f(x)\,dx = \frac{b-a}{2} \int_{-1}^{1} f\left( \frac{a+b}{2} + \frac{b-a}{2}t \right)dt. \]

Całkę na przedziale standardowym przybliża się jako:

\[ \int_{-1}^{1}g(t)\,dt \approx \sum_{i=1}^{m}w_i g(t_i), \]

gdzie \(t_i\) są specjalnie dobranymi węzłami, a \(w_i\) odpowiadającymi im wagami.

Dla dwóch punktów Gaussa:

Dla trzech punktów:

Węzły i wagi można odczytać z tablic albo wyznaczyć jako miejsca zerowe odpowiednich wielomianów ortogonalnych.

Kwadratura Gaussa wykorzystująca \(m\) punktów jest dokładna dla wszystkich wielomianów stopnia nie większego niż:

\[ 2m-1. \]

Jest to bardzo korzystna własność. Dwa odpowiednio dobrane punkty Gaussa zapewniają dokładność dla wielomianów do trzeciego stopnia, a trzy punkty — do stopnia piątego.

Całkowanie adaptacyjne

Równomierny podział przedziału nie zawsze jest efektywny. Funkcja może być prawie liniowa w jednej części przedziału, a w innej gwałtownie zmieniać wartość lub krzywiznę.

W całkowaniu adaptacyjnym przedział dzieli się gęściej tylko tam, gdzie lokalny błąd jest duży. Na spokojnych fragmentach pozostawia się większe podprzedziały.

W adaptacyjnej metodzie Simpsona porównuje się:

• wynik Simpsona na całym przedziale,
• sumę wyników Simpsona na jego dwóch połowach.

Jeżeli \(S\) oznacza wynik na całym przedziale, a \(S_L+S_R\) sumę wyników na połowach, lokalny błąd można oszacować jako:

\[ \frac{|S_L+S_R-S|}{15}. \]

Jeśli błąd przekracza tolerancję, odpowiedni fragment przedziału dzieli się ponownie.

Całkowanie funkcji zadanej tabelarycznie

Całkowanie numeryczne ma szczególne znaczenie, gdy dysponujemy jedynie tabelą wartości:

\[ (x_0,f_0),\ (x_1,f_1),\ldots,\ (x_n,f_n). \]

Metoda trapezów może być stosowana również wtedy, gdy węzły nie są rozmieszczone równomiernie:

\[ T= \sum_{i=0}^{n-1} \frac{x_{i+1}-x_i}{2} \left[ f_i+f_{i+1} \right]. \]

Klasyczna złożona metoda Simpsona wymaga natomiast równomiernego rozstawienia węzłów i parzystej liczby podprzedziałów.

Metoda prostokątów środkowych wymaga znajomości wartości funkcji w środkach podprzedziałów. Jeżeli tabela ich nie zawiera, wartości te trzeba dodatkowo zmierzyć, obliczyć albo przybliżyć za pomocą interpolacji.

Powrót do całki z funkcji \(e^{x^2}\)

Zastosujmy metodę Simpsona do całki:

\[ I= \int_0^1e^{x^2}\,dx. \]

Kolejne wyniki zbliżają się do wartości:

\[ I\approx1{,}4626517459. \]

Przykład ten pokazuje, że brak elementarnej funkcji pierwotnej nie przeszkadza w praktycznym obliczeniu całki oznaczonej.

Zmniejszanie kroku a dokładność

Zmniejszenie kroku zazwyczaj zwiększa dokładność, ale nie należy robić tego bez ograniczeń. Bardzo mały krok powoduje:

• wzrost liczby obliczeń wartości funkcji,
• większy koszt czasowy,
• kumulowanie błędów zaokrągleń,
• możliwość utraty cyfr znaczących.

Wpływ ograniczonej precyzji reprezentacji liczb został opisany w artykule Arytmetyka zmiennoprzecinkowa – cecha, mantysa i zapis binarny.

Oprócz teoretycznych wzorów na błąd można badać zbieżność przez porównywanie wyników uzyskanych dla kroku \(h\) i \(h/2\). Jeżeli błąd metody jest proporcjonalny do \(h^p\), jego praktyczne oszacowanie ma postać:

\[ |E_{h/2}| \approx \frac{|Q_{h/2}-Q_h|}{2^p-1}. \]

Dla trapezów i prostokątów środkowych \(p=2\), natomiast dla metody Simpsona \(p=4\).

Porównanie najważniejszych metod

Którą metodę wybrać?

• Prostokąty lewe lub prawe są łatwe do zrozumienia i implementacji, ale stosunkowo mało dokładne.
• Prostokąty środkowe są proste i zapewniają dokładność drugiego rzędu.
• Trapezy dobrze nadają się do danych tabelarycznych, również przy nierównomiernym rozmieszczeniu węzłów.
• Metoda Simpsona daje bardzo dobre wyniki dla funkcji gładkich i równomiernego podziału.
• Metoda Romberga jest użyteczna, gdy można wielokrotnie zagęszczać podział i zależy nam na dużej dokładności.
• Kwadratury Gaussa są szczególnie efektywne, gdy koszt obliczenia funkcji jest duży i możemy wybierać położenie węzłów.
• Metody adaptacyjne sprawdzają się dla funkcji, których zachowanie silnie zmienia się w różnych częściach przedziału.

Praktyczne zasady całkowania numerycznego

• Sprawdź ciągłość funkcji i położenie ewentualnych osobliwości.
• Dobierz metodę do dostępnej postaci danych.
• Nie zakładaj, że zwiększenie liczby cyfr wyświetlanego wyniku oznacza większą dokładność.
• Porównuj wyniki dla różnych kroków.
• Kontroluj zarówno błąd teoretyczny, jak i praktyczną stabilizację wyniku.
• Przy metodzie Simpsona pamiętaj o parzystej liczbie podprzedziałów.
• Przy Simpsonie 3/8 liczba podprzedziałów musi być podzielna przez trzy.
• W kwadraturach Gaussa prawidłowo przekształć przedział do \([-1,1]\).
• Dziel przedział w miejscach nieciągłości pochodnych lub gwałtownych zmian funkcji.
• Ustal maksymalną liczbę podziałów i minimalny dopuszczalny krok.

Szersze omówienie błędów, ich oszacowań i wpływu pochodnych na wynik znajduje się w artykule Błędy obliczeń numerycznych – rodzaje, źródła i ocena dokładności.

Całkowanie numeryczne nie wymaga znajomości elementarnej funkcji pierwotnej. Wymaga natomiast świadomego wyboru kwadratury, kroku oraz sposobu kontroli błędu.

Podsumowanie

Całkowanie numeryczne pozwala obliczać przybliżone wartości całek oznaczonych, także wtedy, gdy całkowana funkcja nie ma elementarnej funkcji pierwotnej albo jest znana jedynie w postaci danych tabelarycznych.

Metoda prostokątów zastępuje funkcję wartościami stałymi. Jej wariant środkowy jest dokładniejszy od prostokątów lewych i prawych.

Metoda trapezów wykorzystuje średnią wartości funkcji na końcach podprzedziału. Nie należy jej utożsamiać z metodą punktu środkowego, w której obliczamy wartość funkcji dla średniego argumentu.

Metoda Simpsona zastępuje funkcję na parach podprzedziałów parabolami przechodzącymi przez trzy kolejne węzły. Dzięki błędowi rzędu \(h^4\) zapewnia zwykle bardzo dużą dokładność.

Metody Romberga, Gaussa i kwadratury adaptacyjne pozwalają uzyskiwać jeszcze lepsze wyniki przez eliminowanie składników błędu, optymalny dobór węzłów albo lokalne zagęszczanie podziału.

Wzory zawierające supremum pochodnej dają górne ograniczenia błędu. Rzeczywisty błąd może być znacznie mniejszy, co potwierdza przykład całki z funkcji \(1/(1+x^2)\).

Utworzono: 16.06.2026 | Zmodyfikowano: 25.07.2026

Numeryczne rozwiązywanie równań nieliniowych

Mgr inż.

Sebastian Dziarmaga-Działyński

Wirtualna kancelaria korepetytorska i konsultacyjna – wszechwiedza.pl

Wiele problemów matematycznych, technicznych i naukowych prowadzi do równań, których nie można rozwiązać za pomocą prostych przekształceń algebraicznych. Niekiedy rozwiązanie dokładne nie istnieje w postaci elementarnego wzoru, a czasami jego wyznaczenie byłoby zbyt skomplikowane lub niepraktyczne.

W takich sytuacjach stosuje się numeryczne metody rozwiązywania równań. Nie wyznaczają one zazwyczaj pierwiastka w sposób symboliczny, lecz tworzą ciąg kolejnych przybliżeń coraz bliższych poszukiwanemu rozwiązaniu.

W artykule omówimy cztery podstawowe metody rozwiązywania równań nieliniowych:

• metodę bisekcji,
• metodę regula falsi,
• metodę siecznych,
• metodę Newtona.

Różnią się one szybkością zbieżności, wymaganiami dotyczącymi punktów początkowych, koniecznością obliczania pochodnej oraz gwarancją otrzymania rozwiązania.

Podstawowe znaczenie algorytmów przybliżonych, zbieżności i stabilności zostało omówione w artykule Metody numeryczne – czym są i do czego służą?.

Czym jest równanie nieliniowe?

Równanie jednej zmiennej można najczęściej sprowadzić do postaci:

\[f(x)=0.\]

Liczbę \(\alpha\), dla której:

\[f(\alpha)=0,\]

nazywamy pierwiastkiem albo miejscem zerowym funkcji \(f\).

Równanie jest nieliniowe, jeżeli niewiadoma nie występuje wyłącznie w pierwszej potędze i równania nie można zapisać w postaci:

\[ax+b=0.\]

Przykładami równań nieliniowych są:

\[x^3-x-2=0,\]

\[\cos x-x=0,\]

\[e^{-x}-x=0,\]

\[x\sin x-1=0.\]

Rozwiązanie równania \(f(x)=0\) można interpretować geometrycznie jako współrzędną punktu przecięcia wykresu funkcji \(y=f(x)\) z osią \(x\).

Dlaczego potrzebujemy metod numerycznych?

Dla równania liniowego rozwiązanie otrzymujemy bezpośrednio:

\[ax+b=0\quad\Longrightarrow\quad x=-\frac{b}{a}.\]

Istnieją również wzory rozwiązujące równania kwadratowe, sześcienne i równania czwartego stopnia. Nie ma jednak ogólnego wzoru wyrażonego za pomocą działań algebraicznych i pierwiastników dla wszystkich równań wielomianowych stopnia piątego i wyższych.

Jeszcze częściej spotyka się równania zawierające funkcje wykładnicze, logarytmiczne lub trygonometryczne, których pierwiastków nie można przedstawić w prostej postaci dokładnej.

Metody numeryczne pozwalają wtedy znaleźć wartość:

\[\widetilde{\alpha}\approx\alpha\]

z dokładnością wystarczającą do konkretnego zastosowania.

Lokalizacja pierwiastka

Przed rozpoczęciem iteracji warto ustalić, gdzie może znajdować się pierwiastek. Można to zrobić za pomocą:

• wykresu funkcji,
• tabeli wartości,
• analizy znaków funkcji,
• własności monotoniczności,
• oszacowań wynikających z postaci równania,
• wiedzy o analizowanym problemie fizycznym lub technicznym.

Szczególnie ważne jest twierdzenie o wartości pośredniej. Jeżeli funkcja \(f\) jest ciągła na przedziale \([a,b]\) oraz:

\[f(a)\cdot f(b)\lt 0,\]

to w przedziale \((a,b)\) istnieje co najmniej jeden pierwiastek równania:

\[f(x)=0.\]

Warunek zmiany znaku wykorzystują metoda bisekcji i metoda regula falsi.

Warunek ten jest wystarczający, ale nie jest konieczny. Funkcja może dotykać osi \(x\), nie zmieniając znaku. Dzieje się tak na przykład dla pierwiastka parzystej krotności.

Dla funkcji:

\[f(x)=(x-1)^2\]

liczba \(x=1\) jest pierwiastkiem, mimo że po obu jego stronach funkcja pozostaje nieujemna.

Wspólny przykład obliczeniowy

Działanie omawianych metod porównamy na równaniu:

\[x^3-x-2=0.\]

Definiujemy funkcję:

\[f(x)=x^3-x-2.\]

Obliczamy wartości na końcach przedziału \([1,2]\):

\[f(1)=1-1-2=-2,\]

\[f(2)=8-2-2=4.\]

Ponieważ:

\[f(1)\cdot f(2)=-8\lt 0,\]

w przedziale \((1,2)\) znajduje się co najmniej jeden pierwiastek. Jego wartość z większą dokładnością wynosi około:

\[\alpha\approx 1{,}5213797068.\]

Metoda bisekcji

Metoda bisekcji, nazywana również metodą połowienia przedziału, jest jedną z najprostszych i najbardziej niezawodnych metod rozwiązywania równań nieliniowych.

Zakładamy, że funkcja \(f\) jest ciągła na przedziale \([a,b]\) oraz:

\[f(a)\cdot f(b)\lt 0.\]

Obliczamy środek przedziału:

\[c=\frac{a+b}{2}.\]

Następnie sprawdzamy znak \(f(c)\). Zachowujemy tę połowę przedziału, na której końcach funkcja ma przeciwne znaki.

Algorytm bisekcji

1. Wybierz przedział \([a,b]\), dla którego \(f(a)f(b)\lt 0\).
2. Oblicz środek przedziału:
   \[
   c=\frac{a+b}{2}.
   \]
3. Oblicz \(f(c)\).
4. Jeżeli \(f(c)=0\), zakończ obliczenia.
5. Jeżeli \(f(a)f(c)\lt 0\), przyjmij \(b=c\).
6. W przeciwnym przypadku przyjmij \(a=c\).
7. Powtarzaj dzielenie przedziału aż do osiągnięcia wymaganej dokładności.

Przykład bisekcji

Dla równania:

\[x^3-x-2=0\]

rozpoczynamy od przedziału:

\[[a_0,b_0]=[1,2].\]

W każdej iteracji długość przedziału zmniejsza się dwukrotnie. Kolejne przedziały coraz dokładniej ograniczają położenie pierwiastka.

Zbieżność i oszacowanie błędu

Po \(n\) podziałach długość przedziału wynosi:

\[b_n-a_n=\frac{b_0-a_0}{2^n}.\]

Jeżeli jako przybliżenie przyjmiemy środek przedziału:

\[c_n=\frac{a_n+b_n}{2},\]

to błąd spełnia oszacowanie:

\[|\alpha-c_n|\leq\frac{b_0-a_0}{2^{n+1}}.\]

Aby zagwarantować błąd nie większy niż \(\varepsilon\), należy wykonać co najmniej:

\[n\geq \left\lceil \log_2\left(\frac{b_0-a_0}{2\varepsilon}\right) \right\rceil\]

podziałów.

Bisekcja ma zbieżność liniową. Nie jest bardzo szybka, ale pozwala z góry oszacować liczbę potrzebnych iteracji.

Zalety i ograniczenia bisekcji

• Zalety: prosta konstrukcja, gwarantowana zbieżność przy spełnieniu założeń, łatwe oszacowanie błędu i brak konieczności obliczania pochodnej.
• Ograniczenia: stosunkowo wolna zbieżność, konieczność znalezienia przedziału ze zmianą znaku oraz trudność wykrywania pierwiastków parzystej krotności.

Metoda regula falsi

Metoda regula falsi, nazywana metodą fałszywego położenia, podobnie jak bisekcja rozpoczyna się od przedziału \([a,b]\), na którego końcach funkcja przyjmuje wartości przeciwnych znaków.

Zamiast wybierać środek przedziału, prowadzimy prostą przechodzącą przez punkty:

\[(a,f(a))\qquad\text{oraz}\qquad(b,f(b)).\]

Prosta ta jest sieczną wykresu funkcji. Jej punkt przecięcia z osią \(x\) przyjmujemy jako nowe przybliżenie pierwiastka.

Wzór na to przybliżenie ma postać:

\[c= b-\frac{f(b)(b-a)}{f(b)-f(a)}.\]

Równoważnie:

\[c= \frac{af(b)-bf(a)}{f(b)-f(a)}.\]

Algorytm metody regula falsi

1. Wybierz przedział \([a,b]\), dla którego \(f(a)f(b)\lt 0\).
2. Wyznacz punkt przecięcia siecznej z osią \(x\):
   \[
   c=b-\frac{f(b)(b-a)}{f(b)-f(a)}.
   \]
3. Oblicz \(f(c)\).
4. Jeżeli spełnione jest kryterium stopu, zakończ obliczenia.
5. Jeżeli \(f(a)f(c)\lt 0\), przyjmij \(b=c\).
6. W przeciwnym przypadku przyjmij \(a=c\).
7. Powtarzaj obliczenia dla nowego przedziału.

Przykład metody regula falsi

Dla funkcji:

\[f(x)=x^3-x-2\]

oraz przedziału \([1,2]\) pierwsze przybliżenie wynosi:

\[c_1= 2-\frac{4(2-1)}{4-(-2)} = 2-\frac{4}{6} = 1{,}333333.\]

Ponieważ:

\[f(c_1)\approx-0{,}962963,\]

nowy przedział ma postać:

\[[1{,}333333,\ 2].\]

Zbieżność metody regula falsi

Metoda zachowuje przedział zawierający pierwiastek, dzięki czemu przy odpowiednich założeniach jest bezpieczniejsza od metod otwartych. Jej zbieżność jest zazwyczaj liniowa.

Często działa szybciej niż bisekcja, ponieważ wykorzystuje wartości funkcji, a nie tylko ich znaki. Może jednak wystąpić sytuacja, w której jeden koniec przedziału przez wiele iteracji prawie się nie zmienia. Powoduje to spowolnienie obliczeń, nazywane niekiedy zastojem metody.

Zalety i ograniczenia regula falsi

• Zalety: zachowuje przedział zawierający pierwiastek, nie wymaga pochodnej i często jest szybsza od bisekcji.
• Ograniczenia: może zbiegać nierównomiernie, jeden koniec przedziału może pozostawać niemal nieruchomy, a zbieżność nadal jest zazwyczaj liniowa.

Metoda siecznych

Metoda siecznych również wykorzystuje prostą przechodzącą przez dwa punkty wykresu funkcji. W przeciwieństwie do metody regula falsi nie zachowuje jednak przedziału ze zmianą znaku.

Do rozpoczęcia obliczeń potrzebne są dwa przybliżenia:

\[x_0\qquad\text{oraz}\qquad x_1.\]

Sieczna przechodząca przez punkty:

\[(x_{k-1},f(x_{k-1}))\]

oraz:

\[(x_k,f(x_k))\]

przecina oś argumentów w kolejnym przybliżeniu:

\[x_{k+1} = x_k- \frac{f(x_k)(x_k-x_{k-1})} {f(x_k)-f(x_{k-1})}.\]

Wzór można interpretować jako modyfikację metody Newtona, w której pochodna jest zastępowana ilorazem różnicowym:

\[f'(x_k)\approx \frac{f(x_k)-f(x_{k-1})} {x_k-x_{k-1}}.\]

Algorytm metody siecznych

1. Wybierz dwa początkowe przybliżenia \(x_0\) i \(x_1\).
2. Oblicz:
   \[
   x_{k+1}
   =
   x_k-
   \frac{f(x_k)(x_k-x_{k-1})}
   {f(x_k)-f(x_{k-1})}.
   \]
3. Sprawdź kryterium stopu.
4. Przyjmij:
   \[
   x_{k-1}=x_k,\qquad x_k=x_{k+1}.
   \]
5. Powtarzaj obliczenia aż do osiągnięcia wymaganej dokładności.

Przykład metody siecznych

Dla równania:

\[x^3-x-2=0\]

przyjmujemy:

\[x_0=1,\qquad x_1=2.\]

Pierwsze nowe przybliżenie wynosi:

\[x_2 = 2- \frac{4(2-1)} {4-(-2)} = 1{,}333333.\]

W kolejnej iteracji korzystamy już z punktów \(x_1=2\) oraz \(x_2=1{,}333333\), niezależnie od tego, czy tworzą przedział ze zmianą znaku.

Zbieżność metody siecznych

W pobliżu prostego pierwiastka metoda siecznych ma zbieżność superliniową rzędu:

\[\varphi=\frac{1+\sqrt{5}}{2}\approx1{,}618.\]

Jest więc zwykle szybsza od bisekcji i metody regula falsi, ale wolniejsza od prawidłowo działającej metody Newtona.

Nie ma jednak ogólnej gwarancji zbieżności. Przy źle dobranych punktach początkowych kolejne przybliżenia mogą oddalać się od pierwiastka, oscylować albo prowadzić do dzielenia przez bardzo małą wartość:

\[f(x_k)-f(x_{k-1})\approx0.\]

Zalety i ograniczenia metody siecznych

• Zalety: nie wymaga jawnego obliczania pochodnej, ma zbieżność superliniową i jest stosunkowo prosta do zaimplementowania.
• Ograniczenia: wymaga dwóch punktów początkowych, nie gwarantuje zachowania przedziału z pierwiastkiem i może być wrażliwa na wybór wartości początkowych.

Metoda Newtona

Metoda Newtona, nazywana również metodą stycznych, należy do najważniejszych i najszybszych metod rozwiązywania równań nieliniowych.

W punkcie bieżącego przybliżenia \(x_k\) prowadzimy styczną do wykresu funkcji. Jej równanie ma postać:

\[y=f(x_k)+f'(x_k)(x-x_k).\]

Miejsce przecięcia stycznej z osią \(x\) otrzymujemy po przyjęciu \(y=0\):

\[0=f(x_k)+f'(x_k)(x_{k+1}-x_k).\]

Po przekształceniu otrzymujemy wzór Newtona:

\[x_{k+1} = x_k-\frac{f(x_k)}{f'(x_k)}.\]

Metoda wykorzystuje więc liniowe przybliżenie funkcji w otoczeniu aktualnego punktu. Związek stycznej z pochodną został przedstawiony w artykule Pochodna funkcji.

Algorytm metody Newtona

1. Wybierz początkowe przybliżenie \(x_0\).
2. Oblicz wartość funkcji \(f(x_k)\) i pochodnej \(f'(x_k)\).
3. Wyznacz:
   \[
   x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}.
   \]
4. Sprawdź kryterium stopu.
5. Przyjmij \(x_k=x_{k+1}\) i powtórz obliczenia.

Przykład metody Newtona

Dla funkcji:

\[f(x)=x^3-x-2\]

pochodna wynosi:

\[f'(x)=3x^2-1.\]

Wzór iteracyjny ma więc postać:

\[x_{k+1} = x_k- \frac{x_k^3-x_k-2} {3x_k^2-1}.\]

Przyjmijmy wartość początkową:

\[x_0=1{,}5.\]

Pierwsza iteracja daje:

\[x_1 = 1{,}5- \frac{1{,}5^3-1{,}5-2} {3\cdot1{,}5^2-1} \approx1{,}5217391304.\]

Kolejne przybliżenia wynoszą:

W tym przykładzie już po trzech iteracjach otrzymujemy bardzo dokładne przybliżenie pierwiastka.

Zbieżność metody Newtona

Jeżeli punkt początkowy leży dostatecznie blisko prostego pierwiastka \(\alpha\), funkcja jest odpowiednio gładka oraz:

\[f'(\alpha)\ne0,\]

metoda Newtona ma zbieżność kwadratową. Oznacza to w przybliżeniu, że liczba poprawnych cyfr może podwajać się w każdej iteracji.

Jest to jednak wynik lokalny. Metoda nie gwarantuje zbieżności dla dowolnego punktu startowego.

Kiedy metoda Newtona może zawieść?

Problemy mogą wystąpić, gdy:

• punkt początkowy znajduje się daleko od pierwiastka,
• pochodna jest równa lub bardzo bliska zeru,
• funkcja ma wiele pierwiastków i iteracja przechodzi do innego z nich,
• kolejne przybliżenia wpadają w cykl,
• funkcja lub jej pochodna nie jest określona w którymś punkcie iteracji,
• pierwiastek ma krotność większą niż jeden.

Jeżeli:

\[f'(x_k)\approx0,\]

iloraz:

\[\frac{f(x_k)}{f'(x_k)}\]

może przyjąć bardzo dużą wartość i przenieść kolejne przybliżenie daleko od poszukiwanego pierwiastka.

Pierwiastki wielokrotne

Dla pierwiastka wielokrotnego klasyczna metoda Newtona traci zwykle zbieżność kwadratową i staje się zbieżna liniowo.

Jeżeli znana jest krotność \(m\), można zastosować zmodyfikowany wzór:

\[x_{k+1} = x_k- m\frac{f(x_k)}{f'(x_k)}.\]

Zalety i ograniczenia metody Newtona

• Zalety: bardzo szybka zbieżność w pobliżu prostego pierwiastka i wysoka dokładność przy niewielkiej liczbie iteracji.
• Ograniczenia: konieczność obliczania pochodnej, brak globalnej gwarancji zbieżności i duża zależność od wyboru punktu początkowego.

Metody przedziałowe i metody otwarte

Omówione algorytmy można podzielić na dwie grupy.

Metody przedziałowe

• metoda bisekcji,
• metoda regula falsi.

W każdej iteracji zachowują one przedział, na którego końcach funkcja ma przeciwne znaki. Są dzięki temu stosunkowo bezpieczne, ale ich zbieżność jest zazwyczaj liniowa.

Metody otwarte

• metoda siecznych,
• metoda Newtona.

Nie wymagają one zachowywania przedziału obejmującego pierwiastek. Mogą być znacznie szybsze, ale przy niewłaściwym starcie mogą się rozbiegać.

Kryteria zakończenia iteracji

Obliczeń nie prowadzi się w nieskończoność. Iterację kończymy po osiągnięciu zadanej tolerancji albo po przekroczeniu maksymalnej liczby kroków.

Mała wartość funkcji

Można zakończyć obliczenia, gdy:

\[|f(x_k)|\lt\varepsilon_f.\]

Warunek ten sprawdza, czy otrzymany punkt prawie spełnia równanie.

Mała wartość funkcji nie zawsze oznacza jednak mały błąd położenia pierwiastka. Jeżeli wykres jest bardzo płaski, duża zmiana \(x\) może powodować tylko niewielką zmianę \(f(x)\).

Mała zmiana kolejnych przybliżeń

Drugie kryterium ma postać:

\[|x_{k+1}-x_k|\lt\varepsilon_x.\]

Lepsze jest często połączenie tolerancji bezwzględnej i względnej:

\[|x_{k+1}-x_k| \lt \varepsilon_{\mathrm{abs}} + \varepsilon_{\mathrm{rel}}|x_{k+1}|.\]

Szerokość przedziału

W metodach przedziałowych można zastosować warunek:

\[|b_k-a_k|\lt\varepsilon_x.\]

Jest to szczególnie wiarygodne kryterium, ponieważ pierwiastek nadal znajduje się wewnątrz kontrolowanego przedziału.

Maksymalna liczba iteracji

Każdy algorytm powinien mieć również ustaloną maksymalną liczbę iteracji. Chroni to program przed nieskończonym wykonywaniem obliczeń w przypadku rozbieżności, stagnacji lub niewłaściwych danych.

W praktyce dobrze jest stosować jednocześnie kilka kryteriów stopu, na przykład małą wartość funkcji, małą zmianę przybliżeń i ograniczenie liczby iteracji.

Porównanie metod

Którą metodę wybrać?

Nie istnieje jedna metoda najlepsza dla wszystkich równań.

• Bisekcję warto wybrać, gdy najważniejsza jest niezawodność i znany jest przedział ze zmianą znaku.
• Regula falsi może być dobrym wyborem, gdy chcemy zachować bezpieczeństwo metody przedziałowej, ale wykorzystać również wartości funkcji.
• Metoda siecznych jest użyteczna, gdy obliczanie pochodnej jest trudne lub kosztowne, a zależy nam na szybszej zbieżności.
• Metoda Newtona jest zwykle najlepsza, gdy można łatwo obliczyć pochodną i dysponujemy dobrym punktem początkowym.

W praktycznych programach często łączy się różne metody. Najpierw stosuje się bezpieczną bisekcję, aby zawęzić przedział, a następnie przechodzi do szybszej metody Newtona lub siecznych.

Metody hybrydowe

Metoda hybrydowa łączy bezpieczeństwo metody przedziałowej z szybkością metody otwartej.

Przykładowy schemat może wyglądać następująco:

1. Znajdź przedział \([a,b]\), dla którego \(f(a)f(b)\lt0\).
2. Wykonaj kilka kroków bisekcji.
3. Wybierz środek zwężonego przedziału jako punkt startowy metody Newtona.
4. Jeżeli krok Newtona wychodzi poza przedział albo pochodna jest zbyt mała, wykonaj ponownie krok bisekcji.
5. Kontynuuj aż do osiągnięcia zadanej tolerancji.

Takie rozwiązanie ogranicza ryzyko rozbieżności, a jednocześnie pozwala korzystać z szybkiej zbieżności w pobliżu pierwiastka.

Błędy i stabilność obliczeń

Każde przybliżenie pierwiastka jest obarczone błędem. Na wynik wpływają między innymi:

• błąd wynikający z przerwania iteracji,
• błędy zaokrągleń arytmetyki komputerowej,
• niedokładność wartości funkcji i pochodnej,
• uwarunkowanie samego równania,
• niestabilność zastosowanego algorytmu.

Zagadnienia błędu bezwzględnego, względnego, propagacji błędów i uwarunkowania zostały szerzej przedstawione w artykule Błędy obliczeń numerycznych – rodzaje, źródła i ocena dokładności.

W obliczeniach należy także uwzględniać ograniczoną reprezentację liczb opisaną w artykule Arytmetyka zmiennoprzecinkowa – cecha, mantysa i zapis binarny.

Praktyczne zasady implementacji

• Przed rozpoczęciem iteracji sprawdź dziedzinę funkcji.
• W metodach przedziałowych sprawdź ciągłość i zmianę znaku.
• Nie dziel przez wartość pochodnej lub różnicy funkcji bliską zeru.
• Stosuj jednocześnie tolerancję bezwzględną i względną.
• Ustal maksymalną liczbę iteracji.
• Po zakończeniu oblicz wartość reszty \(|f(x_k)|\).
• Sprawdź, czy otrzymane rozwiązanie należy do oczekiwanego zakresu.
• W przypadku wielu pierwiastków wykonuj obliczenia dla różnych przedziałów lub punktów startowych.
• Kontroluj, czy kolejne przybliżenia są liczbami skończonymi.
• Wynik podawaj z liczbą cyfr odpowiadającą rzeczywistej dokładności.

Szybka metoda nie zawsze jest metodą najlepszą. W obliczeniach numerycznych równie ważne jak tempo zbieżności są niezawodność, kontrola błędu i właściwy wybór punktów początkowych.

Podsumowanie

Numeryczne rozwiązywanie równania \(f(x)=0\) polega na konstruowaniu ciągu przybliżeń zmierzających do poszukiwanego pierwiastka.

Metoda bisekcji wielokrotnie dzieli przedział na połowy. Jest wolna, ale niezawodna i pozwala łatwo kontrolować błąd.

Metoda regula falsi zastępuje środek przedziału punktem przecięcia siecznej z osią \(x\), zachowując przy tym przedział obejmujący pierwiastek.

Metoda siecznych korzysta z dwóch ostatnich przybliżeń i nie wymaga obliczania pochodnej. Ma zbieżność superliniową, ale nie gwarantuje otrzymania rozwiązania.

Metoda Newtona wykorzystuje styczną do wykresu funkcji. W pobliżu prostego pierwiastka jest zbieżna kwadratowo, lecz wymaga pochodnej i odpowiedniego punktu startowego.

W praktyce najlepsze rezultaty często dają algorytmy hybrydowe, łączące bezpieczne ograniczanie przedziału z szybką iteracją Newtona lub metodą siecznych.

Utworzono: 15.06.2026 | Zmodyfikowano: 25.07.2026

Błędy obliczeń numerycznych – rodzaje, źródła i ocena dokładności

Mgr inż.

Sebastian Dziarmaga-Działyński

Wirtualna kancelaria korepetytorska i konsultacyjna – wszechwiedza.pl

Wyniki uzyskiwane za pomocą metod numerycznych są najczęściej wynikami przybliżonymi. Różnią się one od wartości dokładnych wskutek niedokładności danych, uproszczeń modelu matematycznego, zastosowania przybliżonego algorytmu oraz ograniczonej precyzji obliczeń komputerowych.

Sam fakt, że wynik jest przybliżony, nie oznacza jednak, że jest on niewiarygodny. Najważniejsze jest określenie wielkości błędu, rozpoznanie jego źródła oraz sprawdzenie, czy otrzymana dokładność jest wystarczająca do danego zastosowania.

Analiza błędów jest jednym z podstawowych zagadnień metod numerycznych. Pozwala ocenić zarówno jakość danych i modelu, jak i zbieżność oraz stabilność zastosowanego algorytmu.

Wartość dokładna i wartość przybliżona

Załóżmy, że dokładna wartość pewnej wielkości wynosi \(x\), natomiast w wyniku pomiaru lub obliczeń otrzymano wartość przybliżoną \(\widetilde{x}\).

Możemy wówczas zapisać:

\[x=\widetilde{x}+e,\]

gdzie \(e\) oznacza błąd ze znakiem:

\[e=x-\widetilde{x}.\]

Jeżeli \(\widetilde{x}x\), błąd ma znak ujemny. W praktyce częściej interesuje nas jednak sama wielkość odchylenia, niezależnie od jego kierunku. Z tego powodu stosuje się błąd bezwzględny.

Błąd bezwzględny

Błąd bezwzględny jest modułem różnicy między wartością dokładną i przybliżoną:

\[\Delta x=\left|x-\widetilde{x}\right|.\]

Błąd bezwzględny ma tę samą jednostkę co analizowana wielkość. Jeżeli na przykład dokładna długość wynosi \(10\ \text{cm}\), a wynik pomiaru jest równy \(9{,}8\ \text{cm}\), to:

\[\Delta x=|10-9{,}8|=0{,}2\ \text{cm}.\]

Wynik można zapisać w postaci:

\[x\approx 9{,}8\ \text{cm},\qquad \Delta x=0{,}2\ \text{cm}.\]

Jeżeli nie znamy wartości dokładnej, często dysponujemy jedynie oszacowaniem maksymalnego błędu. Zapis

\[x=\widetilde{x}\pm\Delta x\]

oznacza wówczas, że poszukiwana wartość powinna należeć do przedziału:

\[x\in[\widetilde{x}-\Delta x,\ \widetilde{x}+\Delta x].\]

Błąd względny i procentowy

Sam błąd bezwzględny nie zawsze pozwala właściwie ocenić dokładność wyniku. Odchylenie równe \(0{,}1\) ma niewielkie znaczenie dla liczby \(1000\), ale jest bardzo duże w porównaniu z liczbą \(0{,}2\).

Błąd względny określa stosunek błędu bezwzględnego do modułu wartości dokładnej:

\[\delta x=\frac{\Delta x}{|x|}=\frac{|x-\widetilde{x}|}{|x|},\qquad x\ne 0.\]

Dla wcześniejszego przykładu:

\[\delta x=\frac{0{,}2}{10}=0{,}02.\]

Błąd procentowy otrzymujemy po pomnożeniu błędu względnego przez \(100\%\):

\[\delta_{\%}x=100\%\cdot\delta x.\]

W rozpatrywanym przypadku:

\[\delta_{\%}x=2\%.\]

Jeżeli wartość dokładna nie jest znana, w mianowniku stosuje się niekiedy wartość przybliżoną:

\[\widetilde{\delta}x=\frac{\Delta x}{|\widetilde{x}|}.\]

Należy wtedy wyraźnie zaznaczyć, że jest to jedynie oszacowanie błędu względnego.

Dokładność i precyzja

Pojęcia dokładności i precyzji bywają używane zamiennie, ale nie oznaczają tego samego.

• Dokładność określa, jak blisko wartości prawdziwej znajduje się wynik.
• Precyzja opisuje szczegółowość zapisu albo powtarzalność wyników.

Seria pomiarów może być bardzo precyzyjna, jeżeli otrzymane wartości są do siebie podobne, ale jednocześnie niedokładna, jeżeli wszystkie są przesunięte względem wartości właściwej.

Podobnie zapisanie wyniku z wieloma cyframi po przecinku nie gwarantuje jego dokładności. Jeżeli dane wejściowe są znane jedynie z dokładnością do dwóch cyfr znaczących, podawanie kilkunastu cyfr wyniku nie zwiększa rzeczywistej ilości informacji.

Główne źródła błędów obliczeń numerycznych

Błąd końcowego wyniku może powstawać na różnych etapach rozwiązywania problemu. Nie wszystkie błędy mają charakter czysto obliczeniowy. Część z nich pojawia się jeszcze przed rozpoczęciem działania algorytmu.

Błędy danych wejściowych

Dane używane w obliczeniach mogą pochodzić z pomiarów, obserwacji, wcześniejszych obliczeń albo tablic. Ich wartości są najczęściej znane tylko z ograniczoną dokładnością.

Źródłem błędu może być między innymi:

• ograniczona dokładność przyrządu pomiarowego,
• wpływ warunków zewnętrznych,
• zaokrąglenie odczytu,
• niewłaściwa kalibracja urządzenia,
• błąd osoby wykonującej pomiar,
• wykorzystanie danych przybliżonych lub nieaktualnych.

Nawet doskonały algorytm nie może na ogół odtworzyć informacji, której nie było w danych wejściowych.

Błędy modelu matematycznego

Model matematyczny jest zawsze uproszczonym opisem rzeczywistości. Podczas jego tworzenia pomija się zjawiska uznane za mniej istotne albo zastępuje skomplikowane zależności prostszymi równaniami.

Przykładowo w modelu ruchu można pominąć opór powietrza, w analizie konstrukcji przyjąć idealną sztywność podpór, a w modelu cieplnym założyć stałość parametrów materiałowych.

Jeżeli model nie opisuje badanego zjawiska wystarczająco dobrze, zwiększanie dokładności obliczeń numerycznych nie usunie błędu modelowania.

Błędy metody

Metoda numeryczna często zastępuje dokładny problem jego prostszym odpowiednikiem. Funkcję można przybliżyć wielomianem, pochodną ilorazem różnicowym, całkę sumą pól, a równanie różniczkowe układem zależności algebraicznych.

Różnica między rozwiązaniem problemu oryginalnego i rozwiązaniem problemu przybliżonego jest nazywana błędem metody lub błędem dyskretyzacji.

Przykładem może być przybliżenie pochodnej wzorem:

\[f'(x)\approx\frac{f(x+h)-f(x)}{h}.\]

Dla skończonego kroku \(h\) prawa strona nie jest na ogół dokładnie równa pochodnej. Błąd maleje zwykle wraz ze zmniejszaniem kroku, ale tylko do momentu, w którym istotne stają się błędy zaokrągleń.

Błędy obcięcia

Błąd obcięcia powstaje wtedy, gdy proces wymagający teoretycznie nieskończonej liczby składników lub iteracji zostaje przerwany po skończonej liczbie kroków.

Przykładowo funkcję wykładniczą można przedstawić za pomocą szeregu:

\[e^x=1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+\ldots\]

W obliczeniach wykorzystuje się jednak jedynie skończoną liczbę początkowych składników:

\[e^x\approx 1+x+\frac{x^2}{2!}+\ldots+\frac{x^n}{n!}.\]

Odrzucona część szeregu stanowi źródło błędu obcięcia.

Błędy zaokrągleń

Komputer zapisuje liczby za pomocą skończonej liczby bitów. Nie wszystkie liczby rzeczywiste można więc przedstawić dokładnie. Wiele wartości zostaje zastąpionych najbliższą liczbą możliwą do zapisania.

Błędy te są szczegółowo związane z arytmetyką zmiennoprzecinkową. Przykładem jest zapis liczb \(0{,}1\) i \(0{,}2\), które mają nieskończone rozwinięcia w systemie binarnym.

Każde kolejne działanie może wprowadzać nowy błąd zaokrąglenia. Przy dużej liczbie operacji błędy mogą się kumulować lub wzajemnie częściowo znosić.

Propagacja błędów

Wielkości wykorzystywane w obliczeniach są często obarczone błędami. Jeżeli wynik zależy od jednej lub wielu takich wielkości, ich niedokładności wpływają na niedokładność wyniku końcowego. Zjawisko to nazywa się propagacją błędów.

Dodawanie i odejmowanie

Niech:

\[z=x\pm y.\]

Jeżeli maksymalne błędy bezwzględne argumentów wynoszą \(\Delta x\) i \(\Delta y\), to maksymalny błąd bezwzględny wyniku można oszacować jako:

\[\Delta z\leq\Delta x+\Delta y.\]

Przy odejmowaniu błędy bezwzględne również się sumują. Szczególnie niebezpieczne jest odejmowanie liczb prawie równych, ponieważ mały wynik może mieć bardzo duży błąd względny.

Mnożenie i dzielenie

Dla iloczynu:

\[z=xy\]

w przybliżeniu pierwszego rzędu otrzymujemy:

\[\delta z\approx\delta x+\delta y.\]

Podobna zależność obowiązuje dla ilorazu:

\[z=\frac{x}{y},\qquad \delta z\approx\delta x+\delta y.\]

W przypadku szacowania najgorszego możliwego błędu względnego dodaje się wartości bezwzględne poszczególnych wkładów.

Potęgowanie

Jeżeli:

\[z=x^\alpha,\]

to dla niewielkich błędów:

\[\delta z\approx|\alpha|\delta x.\]

Przykładowo względny błąd kwadratu jest w przybliżeniu dwukrotnie większy od względnego błędu podstawy:

\[\delta(x^2)\approx 2\delta x.\]

Propagacja błędów metodą różniczki zupełnej

Ogólny sposób szacowania wpływu niewielkich błędów danych na wynik opiera się na rachunku różniczkowym. Jeżeli:

\[y=f(x_1,x_2,\ldots,x_n),\]

to różniczka zupełna funkcji ma postać:

\[dy=\frac{\partial f}{\partial x_1}dx_1+\frac{\partial f}{\partial x_2}dx_2+\ldots+\frac{\partial f}{\partial x_n}dx_n.\]

Pochodna cząstkowa

\[\frac{\partial f}{\partial x_i}\]

określa wrażliwość wyniku na zmianę zmiennej \(x_i\), gdy pozostałe zmienne są traktowane jako stałe. Podstawowe znaczenie pochodnej zostało omówione w artykule Pochodna funkcji.

Maksymalny błąd bezwzględny

Jeżeli błędy argumentów są niewielkie, można zastąpić różniczki oszacowaniami błędów. W najgorszym możliwym przypadku przyjmujemy sumę modułów wszystkich składników:

\[\Delta y\approx \left|\frac{\partial f}{\partial x_1}\right|\Delta x_1+ \left|\frac{\partial f}{\partial x_2}\right|\Delta x_2+ \ldots+ \left|\frac{\partial f}{\partial x_n}\right|\Delta x_n.\]

Pochodne cząstkowe oblicza się dla przyjętych wartości argumentów. Moduły zapewniają, że poszczególne wkłady nie znoszą się podczas wyznaczania maksymalnego możliwego błędu.

Błąd względny funkcji wielu zmiennych

Po podzieleniu oszacowania przez moduł wartości funkcji otrzymujemy przybliżony błąd względny:

\[\delta y\approx \frac{1}{|f|} \sum_{i=1}^{n} \left|\frac{\partial f}{\partial x_i}\right| \Delta x_i.\]

Jeżeli błędy argumentów zapisujemy jako błędy względne:

\[\Delta x_i=|x_i|\delta x_i,\]

otrzymujemy:

\[\delta y\approx \sum_{i=1}^{n} \left| \frac{x_i}{f} \frac{\partial f}{\partial x_i} \right|\delta x_i.\]

Współczynniki

\[\left| \frac{x_i}{f} \frac{\partial f}{\partial x_i} \right|\]

pokazują, w jakim stopniu względny błąd danej zmiennej wpływa na względny błąd wyniku.

Iloczyny, ilorazy i potęgi

Dla funkcji:

\[y=Cx_1^{\alpha_1}x_2^{\alpha_2}\ldots x_n^{\alpha_n}\]

wzór na maksymalny błąd względny upraszcza się do postaci:

\[\delta y\approx |\alpha_1|\delta x_1+ |\alpha_2|\delta x_2+ \ldots+ |\alpha_n|\delta x_n.\]

Przykładowo dla:

\[y=\frac{x^2z}{w^3}\]

otrzymujemy:

\[\delta y\approx 2\delta x+\delta z+3\delta w.\]

Przykład propagacji błędu – pole prostokąta

Załóżmy, że zmierzono boki prostokąta:

\[a=10{,}0\pm 0{,}1\ \text{cm},\]

\[b=5{,}0\pm 0{,}05\ \text{cm}.\]

Pole wynosi:

\[P=ab.\]

Dla wartości przybliżonych otrzymujemy:

\[P=10{,}0\cdot 5{,}0=50{,}0\ \text{cm}^2.\]

Pochodne cząstkowe są równe:

\[\frac{\partial P}{\partial a}=b,\qquad \frac{\partial P}{\partial b}=a.\]

Maksymalny błąd bezwzględny pola szacujemy jako:

\[\Delta P\approx \left|\frac{\partial P}{\partial a}\right|\Delta a+ \left|\frac{\partial P}{\partial b}\right|\Delta b.\]

Po podstawieniu:

\[\Delta P\approx 5{,}0\cdot 0{,}1+ 10{,}0\cdot 0{,}05 =1{,}0\ \text{cm}^2.\]

Wynik można zapisać jako:

\[P=(50{,}0\pm 1{,}0)\ \text{cm}^2.\]

Maksymalny błąd względny wynosi:

\[\delta P=\frac{1{,}0}{50{,}0}=0{,}02=2\%.\]

Propagacja niezależnych niepewności

Suma modułów opisuje najgorszy możliwy przypadek, w którym wszystkie błędy działają jednocześnie w kierunku zwiększającym błąd wyniku. Jeżeli niedokładności wejściowe mają charakter losowy i są od siebie niezależne, zwykle stosuje się pierwiastek z sumy kwadratów.

Dla funkcji:

\[y=f(x_1,x_2,\ldots,x_n)\]

niepewność standardową wyniku można oszacować ze wzoru:

\[u(y)= \sqrt{ \left(\frac{\partial f}{\partial x_1}u(x_1)\right)^2+ \left(\frac{\partial f}{\partial x_2}u(x_2)\right)^2+ \ldots+ \left(\frac{\partial f}{\partial x_n}u(x_n)\right)^2 }.\]

Dla pola prostokąta:

\[u(P)= \sqrt{ \bigl(bu(a)\bigr)^2+ \bigl(au(b)\bigr)^2 }.\]

Jeżeli przyjmiemy:

\[u(a)=0{,}1\ \text{cm},\qquad u(b)=0{,}05\ \text{cm},\]

to:

\[u(P)= \sqrt{ (5\cdot 0{,}1)^2+ (10\cdot 0{,}05)^2 } \approx 0{,}71\ \text{cm}^2.\]

Wartość ta jest mniejsza od maksymalnego błędu \(1{,}0\ \text{cm}^2\), ponieważ nie zakładamy jednoczesnego wystąpienia wszystkich odchyleń w najbardziej niekorzystnym kierunku.

Wielkości skorelowane

Jeżeli zmienne wejściowe nie są niezależne, należy uwzględnić ich kowariancje. Pełny wzór ma postać:

\[u^2(y)= \sum_{i=1}^{n} \left(\frac{\partial f}{\partial x_i}\right)^2u^2(x_i) + 2\sum_{i<j} \frac{\partial f}{\partial x_i} \frac{\partial f}{\partial x_j} \operatorname{cov}(x_i,x_j).\]

Dodatnia korelacja może zwiększać niepewność wyniku, natomiast korelacja ujemna może powodować częściowe znoszenie wpływu błędów.

Kiedy metoda różniczkowa jest wiarygodna?

Propagacja błędów za pomocą różniczki zupełnej jest przybliżeniem liniowym. Daje dobre wyniki przede wszystkim wtedy, gdy:

• błędy argumentów są niewielkie,
• funkcja jest różniczkowalna w rozpatrywanym otoczeniu,
• pochodne nie zmieniają się gwałtownie,
• punkt obliczeń nie leży w pobliżu osobliwości,
• składniki wyższych rzędów można pominąć.

Jeżeli niepewności są duże albo funkcja jest silnie nieliniowa, dokładniejsze wyniki można otrzymać przez obliczenie wartości skrajnych, rozwinięcie wyższego rzędu albo symulację Monte Carlo.

Uwarunkowanie zadania

Uwarunkowanie opisuje wrażliwość rozwiązania na niewielkie zmiany danych wejściowych. Jest właściwością samego problemu matematycznego, a nie zastosowanego algorytmu.

Dla funkcji jednej zmiennej:

\[y=f(x)\]

względną liczbę uwarunkowania można zapisać jako:

\[\kappa(x)= \left| \frac{x f'(x)}{f(x)} \right|.\]

W przybliżeniu:

\[\delta y\approx\kappa(x)\delta x.\]

Jeżeli \(\kappa(x)\) jest niewielka, małe względne zmiany danych powodują małe względne zmiany wyniku. Jeżeli liczba uwarunkowania jest bardzo duża, zadanie jest źle uwarunkowane.

Źle uwarunkowanego problemu nie można naprawić wyłącznie przez zastosowanie lepszego algorytmu. Można jednak próbować zmienić postać zadania, skalować dane albo zwiększyć precyzję obliczeń.

Stabilność algorytmu

Stabilność jest cechą konkretnego sposobu wykonywania obliczeń. Algorytm stabilny nie wzmacnia nadmiernie błędów zaokrągleń i niewielkich zaburzeń danych.

Należy więc odróżnić:

• źle uwarunkowane zadanie — wynik jest bardzo wrażliwy na dane niezależnie od metody,
• niestabilny algorytm — zastosowany sposób obliczeń niepotrzebnie wzmacnia błędy.

Dla dobrze uwarunkowanego problemu można zwykle znaleźć algorytm dający wiarygodne wyniki. Zastosowanie niestabilnej metody może jednak zniszczyć dokładność nawet wtedy, gdy samo zadanie nie jest szczególnie wrażliwe.

Utrata cyfr znaczących

Jednym z najważniejszych problemów numerycznych jest odejmowanie dwóch prawie równych liczb. Wspólne początkowe cyfry redukują się, a w wyniku pozostają cyfry, które mogły być wcześniej obarczone błędem zaokrąglenia.

Zjawisko to nazywa się utratą cyfr znaczących albo katastrofalnym odejmowaniem.

Rozważmy wyrażenie:

\[\sqrt{x+1}-\sqrt{x}.\]

Dla dużego \(x\) obie wartości podlegające odejmowaniu są do siebie bardzo podobne. Wyrażenie można jednak przekształcić:

\[\sqrt{x+1}-\sqrt{x} = \frac{(\sqrt{x+1}-\sqrt{x})(\sqrt{x+1}+\sqrt{x})} {\sqrt{x+1}+\sqrt{x}},\]

czyli:

\[\sqrt{x+1}-\sqrt{x} = \frac{1}{\sqrt{x+1}+\sqrt{x}}.\]

Druga postać jest matematycznie równoważna, ale zwykle znacznie stabilniejsza numerycznie.

Kumulowanie się błędów

Podczas wykonywania wielu działań błędy mogą się kumulować. Szczególne znaczenie ma kolejność obliczeń.

W arytmetyce liczb rzeczywistych dodawanie jest łączne:

\[(a+b)+c=a+(b+c).\]

W arytmetyce zmiennoprzecinkowej po każdym działaniu następuje zaokrąglenie, dlatego może wystąpić:

\[\operatorname{fl}(\operatorname{fl}(a+b)+c) \ne \operatorname{fl}(a+\operatorname{fl}(b+c)).\]

Przy sumowaniu wielu wartości można poprawić dokładność przez:

• sumowanie od najmniejszego do największego modułu,
• sumowanie parami,
• zastosowanie algorytmu Kahana,
• wykorzystanie większej precyzji,
• odpowiednie skalowanie danych.

Błąd metody a błąd zaokrągleń

Zmniejszanie kroku metody numerycznej zwykle ogranicza błąd dyskretyzacji. Nie oznacza to jednak, że krok należy zmniejszać bez końca.

Bardzo mały krok może powodować:

• wykonywanie znacznie większej liczby działań,
• kumulowanie błędów zaokrągleń,
• odejmowanie bardzo podobnych wartości,
• utratę cyfr znaczących.

Łączny błąd można schematycznie zapisać jako:

\[\text{błąd całkowity} \approx \text{błąd metody} + \text{błąd zaokrągleń}.\]

W wielu problemach istnieje więc optymalny zakres kroku, dla którego łączny błąd jest najmniejszy.

Cyfry znaczące i zapis wyniku

Wynik powinien być zapisywany w sposób zgodny z jego rzeczywistą dokładnością. Jeżeli otrzymano:

\[x=12{,}347826\pm 0{,}2,\]

podawanie wszystkich cyfr wartości przybliżonej jest nieuzasadnione. Właściwszy zapis to:

\[x=12{,}3\pm 0{,}2.\]

Zwykle niepewność zaokrągla się do jednej lub dwóch cyfr znaczących, a wartość przybliżoną do tego samego miejsca dziesiętnego.

Należy również odróżnić liczbę cyfr wyświetlanych przez kalkulator lub program od liczby cyfr, które są rzeczywiście wiarygodne.

Jak ograniczać błędy numeryczne?

• Analizować jakość i dokładność danych wejściowych.
• Sprawdzać, czy model matematyczny jest odpowiedni do badanego zjawiska.
• Dobierać metodę o znanych własnościach zbieżności.
• Stosować algorytmy stabilne numerycznie.
• Unikać odejmowania liczb prawie równych.
• Skalować zmienne o bardzo różnych rzędach wielkości.
• Nie stosować kroku obliczeniowego mniejszego, niż wymaga tego zadanie.
• Porównywać wyniki uzyskane dla różnych kroków i tolerancji.
• Kontrolować błąd bezwzględny i względny.
• W razie potrzeby zwiększać precyzję obliczeń.
• Sprawdzać wynik za pomocą innej metody albo rozwiązania wzorcowego.
• Analizować jednostki i rząd wielkości odpowiedzi.

Wynik numeryczny bez informacji o jego dokładności może wyglądać wiarygodnie, nawet jeżeli jest obarczony dużym błędem.

Błędy w aproksymacji i interpolacji

Analiza błędów jest istotna również podczas przybliżania danych za pomocą funkcji. W interpolacji należy uwzględnić błąd wynikający z wyboru węzłów i stopnia wielomianu, natomiast w aproksymacji ocenia się odchylenia funkcji przybliżającej od obserwowanych wartości.

Zagadnienia te zostały szerzej opisane w artykule Aproksymacja i interpolacja.

Podsumowanie

Błędy obliczeń numerycznych mogą wynikać z niedokładności danych wejściowych, uproszczeń modelu, zastosowanej metody, przerwania nieskończonego procesu oraz ograniczonej precyzji arytmetyki komputerowej.

Błąd bezwzględny określa odległość między wartością dokładną i przybliżoną, natomiast błąd względny odnosi tę różnicę do skali analizowanej wielkości.

Wpływ błędów argumentów na wynik funkcji wielu zmiennych można oszacować za pomocą różniczki zupełnej i pochodnych cząstkowych. Suma modułów opisuje maksymalny możliwy błąd, natomiast pierwiastek z sumy kwadratów jest stosowany dla niezależnych niepewności losowych.

Uwarunkowanie jest właściwością problemu matematycznego, a stabilność — właściwością algorytmu. Rozróżnienie tych pojęć jest niezbędne do właściwej oceny jakości wyników numerycznych.

Świadoma kontrola błędów pozwala dobierać odpowiednią metodę, precyzję i sposób zapisu wyniku oraz uniknąć sytuacji, w której długa liczba wyświetlana przez komputer jedynie pozornie oznacza dużą dokładność.

Utworzono: 14.06.2026 | Zmodyfikowano: 25.07.2026

Arytmetyka zmiennoprzecinkowa – cecha, mantysa i zapis binarny

Mgr inż.

Sebastian Dziarmaga-Działyński

Wirtualna kancelaria korepetytorska i konsultacyjna – wszechwiedza.pl

Komputer wykonuje obliczenia na liczbach zapisanych za pomocą skończonej liczby bitów. Oznacza to, że nie może dokładnie przechowywać wszystkich liczb rzeczywistych. Wiele wartości musi zostać zastąpionych najbliższą liczbą możliwą do przedstawienia w pamięci urządzenia.

Arytmetyka zmiennoprzecinkowa określa sposób zapisywania i przetwarzania takich przybliżonych wartości. Umożliwia reprezentowanie zarówno bardzo dużych, jak i bardzo małych liczb, ale wprowadza błędy zaokrągleń, ograniczoną liczbę cyfr znaczących oraz szczególne przypadki, takie jak nieskończoność i wartość NaN.

Zrozumienie arytmetyki zmiennoprzecinkowej jest szczególnie ważne w metodach numerycznych. Nawet prawidłowo dobrany algorytm może zwrócić niedokładny wynik, jeżeli występujące podczas obliczeń błędy zaokrągleń są wzmacniane przez kolejne działania.

Co oznacza zapis zmiennoprzecinkowy?

Nazwa „zmiennoprzecinkowy” oznacza, że położenie przecinka nie jest ustalone. Liczba jest przechowywana w sposób przypominający zapis wykładniczy stosowany w matematyce.

W systemie dziesiętnym możemy zapisać na przykład:

\[123\,450=1{,}2345\cdot 10^5\]

oraz

\[0{,}0012345=1{,}2345\cdot 10^{-3}.\]

W obu przypadkach część znacząca liczby jest taka sama, natomiast zmienia się wykładnik potęgi dziesięciu. Dzięki temu można zapisywać liczby o bardzo różnym rzędzie wielkości.

Komputer postępuje podobnie, lecz zamiast systemu dziesiętnego wykorzystuje system binarny. Ogólną postać znormalizowanej liczby zmiennoprzecinkowej można zapisać jako

\[x=(-1)^s\cdot m\cdot 2^e,\]

gdzie:

• \(s\) określa znak liczby,
• \(m\) jest częścią znaczącą, nazywaną tradycyjnie mantysą,
• \(e\) jest wykładnikiem potęgi liczby \(2\).

W literaturze informatycznej zamiast słowa „mantysa” często używa się dokładniejszego określenia znacząca lub angielskiego terminu significand. W potocznym opisie formatu IEEE 754 nadal powszechnie mówi się jednak o polu mantysy.

Znormalizowany zapis binarny

Podobnie jak w dziesiętnym zapisie naukowym, liczba binarna może zostać przedstawiona w postaci znormalizowanej. Przecinek przesuwa się tak, aby przed nim znalazła się dokładnie jedna niezerowa cyfra.

W systemie binarnym jedyną niezerową cyfrą jest \(1\), dlatego dodatnia liczba znormalizowana ma postać

\[1{,}b_1b_2b_3\ldots{}_2\cdot 2^e.\]

Przykładowo liczba binarna

\[101011{,}01_2\]

może zostać zapisana jako

\[1{,}0101101_2\cdot 2^5.\]

Przecinek przesunięto o pięć miejsc w lewo, dlatego wykładnik wynosi \(5\). Pierwsza jedynka w znormalizowanym zapisie jest zawsze obecna, więc nie trzeba zapisywać jej w pamięci. Nazywa się ją ukrytą jedynką.

Dzięki ukrytej jedynce liczba przechowywanych bitów części znaczącej jest o jeden mniejsza od rzeczywistej precyzji zapisu. W formacie pojedynczej precyzji przechowuje się 23 bity części ułamkowej, lecz znacząca ma efektywnie 24 bity.

Standard IEEE 754

Najczęściej stosowane zasady reprezentacji liczb zmiennoprzecinkowych określa standard IEEE 754. Definiuje on format zapisu liczb, wykonywanie działań arytmetycznych, sposoby zaokrąglania oraz wartości specjalne.

W typowym zapisie binarnym liczba jest dzielona na trzy pola:

• znak — określa, czy liczba jest dodatnia, czy ujemna,
• cecha — przechowuje wykładnik zapisany z odpowiednim przesunięciem,
• mantysa — przechowuje bity części znaczącej znajdujące się po ukrytej jedynce.

Pojedyncza precyzja

Format pojedynczej precyzji, określany również jako binary32 lub często jako typ float, wykorzystuje 32 bity:

• 1 bit znaku,
• 8 bitów cechy,
• 23 bity mantysy.

Wartość przesunięcia wykładnika wynosi w tym formacie

\[\operatorname{bias}=127.\]

Jeżeli w polu cechy zapisano wartość \(E\), rzeczywisty wykładnik liczby znormalizowanej wynosi

\[e=E-127.\]

Format pojedynczej precyzji zapewnia około 7 cyfr dziesiętnych dokładności. Największa skończona dodatnia liczba ma wartość około

\[3{,}4028235\cdot 10^{38}.\]

Podwójna precyzja

Format podwójnej precyzji, określany jako binary64 lub zwykle jako typ double, wykorzystuje 64 bity:

• 1 bit znaku,
• 11 bitów cechy,
• 52 bity mantysy.

Przesunięcie wykładnika wynosi

\[\operatorname{bias}=1023.\]

Podwójna precyzja zapewnia około 15 – 16 poprawnych cyfr dziesiętnych. Jest najczęściej stosowana w obliczeniach naukowych i numerycznych, ponieważ znacznie ogranicza błędy zaokrągleń w porównaniu z pojedynczą precyzją.

Przykład zapisu liczby 43,25

Rozważmy liczbę dziesiętną

\[43{,}25_{10}.\]

Jej część całkowita ma w systemie binarnym postać

\[43_{10}=101011_2,\]

natomiast część ułamkowa

\[0{,}25_{10}=0{,}01_2.\]

Otrzymujemy zatem

\[43{,}25_{10}=101011{,}01_2.\]

Po normalizacji:

\[101011{,}01_2=1{,}0101101_2\cdot 2^5.\]

Liczba jest dodatnia, dlatego bit znaku wynosi

\[s=0.\]

Rzeczywisty wykładnik jest równy \(5\). W formacie pojedynczej precyzji trzeba dodać do niego przesunięcie 127:

\[E=5+127=132.\]

Liczba \(132\) ma zapis binarny

\[132_{10}=10000100_2.\]

W polu mantysy umieszcza się bity znajdujące się po pierwszej jedynce znormalizowanego zapisu:

\[01011010000000000000000.\]

Pełny 32-bitowy zapis ma więc postać:

\[\underbrace{0}_{\text{znak}}\ \underbrace{10000100}_{\text{cecha}}\ \underbrace{01011010000000000000000}_{\text{mantysa}}.\]

Dlaczego nie wszystkie liczby można zapisać dokładnie?

W systemie dziesiętnym niektóre ułamki mają skończone rozwinięcie, a inne rozwinięcie nieskończone. Przykładowo:

\[\frac{1}{4}=0{,}25,\]

natomiast

\[\frac{1}{3}=0{,}333333\ldots\]

W systemie binarnym ułamek ma skończone rozwinięcie tylko wtedy, gdy po skróceniu jego mianownik jest potęgą liczby \(2\). Dlatego liczby takie jak \(0{,}5\), \(0{,}25\) i \(0{,}125\) można zapisać dokładnie:

\[0{,}5_{10}=0{,}1_2,\]

\[0{,}25_{10}=0{,}01_2,\]

\[0{,}125_{10}=0{,}001_2.\]

Liczby \(0{,}1\) i \(0{,}2\) mają natomiast nieskończone okresowe rozwinięcia binarne:

\[0{,}1_{10}=0{,}00011001100110011\ldots{}_2,\]

\[0{,}2_{10}=0{,}0011001100110011\ldots{}_2.\]

Komputer dysponuje skończoną liczbą bitów, dlatego musi obciąć lub zaokrąglić nieskończone rozwinięcie. W pamięci zostaje zapisana liczba bardzo bliska wartości właściwej, lecz nie zawsze dokładnie jej równa.

Dlaczego \(0{,}1+0{,}2\) nie zawsze daje dokładnie \(0{,}3\)?

Zarówno \(0{,}1\), jak i \(0{,}2\) są przechowywane w komputerze jako przybliżenia. Po wykonaniu dodawania powstaje kolejna liczba przybliżona. W podwójnej precyzji jej dokładna wartość może zostać wyświetlona jako

\[0{,}30000000000000004.\]

Różnica względem \(0{,}3\) jest niezwykle mała, ale może być widoczna po wyświetleniu odpowiednio dużej liczby cyfr. Nie oznacza to uszkodzenia procesora ani błędu języka programowania. Jest to naturalna konsekwencja binarnej reprezentacji liczb.

Zaokrąglanie liczb zmiennoprzecinkowych

Gdy dokładny wynik działania nie może zostać zapisany w wybranym formacie, musi zostać zastąpiony jedną z sąsiednich liczb reprezentowalnych. Standard IEEE 754 definiuje kilka sposobów zaokrąglania.

• Zaokrąglanie do najbliższej liczby — wybierana jest najbliższa wartość reprezentowalna.
• Zaokrąglanie w kierunku zera — część wykraczająca poza dostępną precyzję jest odrzucana.
• Zaokrąglanie w kierunku \(+\infty\) — wybierana jest najmniejsza liczba reprezentowalna nie mniejsza od wyniku dokładnego.
• Zaokrąglanie w kierunku \(-\infty\) — wybierana jest największa liczba reprezentowalna nie większa od wyniku dokładnego.

Domyślnie stosuje się zazwyczaj zaokrąglanie do najbliższej wartości z regułą „do parzystej” w przypadku wyniku leżącego dokładnie w połowie między dwiema liczbami. Ogranicza to systematyczne przesuwanie wyników w jednym kierunku podczas długich serii obliczeń.

Epsilon maszynowy

Epsilon maszynowy określa odległość między liczbą \(1\) a następną większą liczbą możliwą do zapisania w danym formacie zmiennoprzecinkowym.

Można go opisać warunkiem

\[1+\varepsilon_{\mathrm{mach}}>1,\]

przy czym dla odpowiednio mniejszych liczb dodanie ich do jedynki może nie zmienić zapisanego wyniku:

\[1+\frac{\varepsilon_{\mathrm{mach}}}{2}=1.\]

Dla pojedynczej precyzji:

\[\varepsilon_{\mathrm{mach}}=2^{-23}\approx 1{,}1920929\cdot 10^{-7}.\]

Dla podwójnej precyzji:

\[\varepsilon_{\mathrm{mach}}=2^{-52}\approx 2{,}220446049250313\cdot 10^{-16}.\]

W niektórych podręcznikach przez epsilon maszynowy rozumie się maksymalny względny błąd pojedynczego zaokrąglenia. Wówczas stosowana wartość jest dwukrotnie mniejsza i bywa nazywana zaokrągleniem jednostkowym:

\[u=\frac{\varepsilon_{\mathrm{mach}}}{2}.\]

Epsilon maszynowy nie oznacza najmniejszej dodatniej liczby, jaką można zapisać. Opisuje on rozdzielczość liczb znajdujących się w pobliżu jedynki. Odległość między sąsiednimi wartościami rośnie wraz ze wzrostem rzędu wielkości liczby.

Rozmieszczenie liczb zmiennoprzecinkowych

Liczby zmiennoprzecinkowe nie są rozmieszczone na osi liczbowej równomiernie. W pobliżu zera leżą bardzo gęsto, natomiast dla dużych wartości odległość między kolejnymi liczbami rośnie.

W określonym przedziale wykładnika odstęp między sąsiednimi wartościami jest stały. Po przejściu do kolejnej potęgi liczby \(2\) odstęp ten podwaja się.

Oznacza to, że duża liczba może nie zmienić się po dodaniu do niej stosunkowo małej wartości. Przykładowo w ograniczonej precyzji może wystąpić sytuacja

\[x+1=x,\]

jeżeli odległość między sąsiednimi liczbami reprezentowalnymi w pobliżu \(x\) jest większa od \(1\).

Liczby znormalizowane i zdenormalizowane

Większość liczb w formacie IEEE 754 jest zapisywana w postaci znormalizowanej, z ukrytą jedynką przed przecinkiem binarnym:

\[x=(-1)^s\cdot (1{,}f)_2\cdot 2^{E-\operatorname{bias}}.\]

Dla bardzo małych wartości wykładnik nie może być już dalej zmniejszany. Stosuje się wtedy liczby zdenormalizowane, nazywane również subnormalnymi.

Dla liczby subnormalnej ukryta jedynka zostaje zastąpiona zerem:

\[x=(-1)^s\cdot (0{,}f)_2\cdot 2^{1-\operatorname{bias}}.\]

Liczby subnormalne umożliwiają łagodne przejście między najmniejszymi liczbami znormalizowanymi a zerem. Precyzja względna stopniowo maleje, ale bardzo mały wynik nie zostaje natychmiast zastąpiony zerem.

Wartości specjalne

Nie wszystkie kombinacje bitów w standardzie IEEE 754 oznaczają zwykłe liczby skończone. Część z nich przeznaczono do reprezentowania zera, nieskończoności, liczb subnormalnych i wartości nieokreślonych.

Dodatnie i ujemne zero

Standard rozróżnia \(+0\) i \(-0\). W większości porównań wartości te są sobie równe, ale znak zera może mieć znaczenie w niektórych działaniach, na przykład:

\[\frac{1}{+0}=+\infty,\qquad \frac{1}{-0}=-\infty.\]

Nieskończoność

Nieskończoność może powstać między innymi po przekroczeniu największej reprezentowalnej wartości albo podczas dzielenia liczby niezerowej przez zero.

Przykładowo:

\[\frac{1}{0}=+\infty,\]

o ile środowisko obliczeniowe stosuje reguły IEEE 754 i nie przerywa działania wyjątkiem.

Wartość NaN

Skrót NaN pochodzi od angielskiego określenia Not a Number. Wartość ta oznacza wynik operacji, której nie można przypisać zwykłej wartości liczbowej.

NaN może powstać na przykład podczas obliczania:

• \(\frac{0}{0}\),
• \(\infty-\infty\),
• \(\sqrt{-1}\) w arytmetyce liczb rzeczywistych,
• logarytmu liczby ujemnej w arytmetyce liczb rzeczywistych.

Wartość NaN ma nietypowe własności porównawcze. Nie jest równa żadnej liczbie, a nawet samej sobie:

\[\operatorname{NaN}\ne\operatorname{NaN}.\]

Nadmiar i niedomiar

Nadmiar

Nadmiar, nazywany po angielsku overflow, występuje wtedy, gdy wartość wyniku przekracza zakres danego formatu.

W zależności od sposobu wykonywania obliczeń wynikiem może być nieskończoność, największa liczba skończona albo zgłoszenie błędu.

Niedomiar

Niedomiar, czyli underflow, występuje wtedy, gdy dodatnia wartość jest zbyt mała, aby można ją było przedstawić jako liczbę znormalizowaną.

Początkowo wynik może zostać zapisany jako liczba subnormalna. Po dalszym zmniejszaniu zostaje zaokrąglony do zera.

Najczęstsze problemy w obliczeniach

Działania nie zawsze są łączne

W matematyce rzeczywistej dodawanie jest działaniem łącznym:

\[(a+b)+c=a+(b+c).\]

W arytmetyce zmiennoprzecinkowej równość ta nie musi być spełniona, ponieważ po każdym działaniu następuje zaokrąglenie.

Jeżeli \(a\) jest bardzo duże, a \(b\) i \(c\) znacznie mniejsze, może wystąpić sytuacja:

\[(a+b)+c\ne a+(b+c).\]

Kolejność sumowania może więc wpływać na wynik. Przy dodawaniu wielu liczb często korzystniej jest rozpoczynać od wartości o najmniejszym module.

Utrata cyfr znaczących

Szczególnie niebezpieczne jest odejmowanie dwóch prawie równych liczb. Ich początkowe cyfry znaczące redukują się, a w wyniku pozostają głównie cyfry obarczone błędem zaokrąglenia.

Zjawisko to nazywa się katastrofalnym odejmowaniem lub utratą cyfr znaczących.

Przykładowo wyrażenie

\[\sqrt{x+1}-\sqrt{x}\]

dla dużego \(x\) prowadzi do odejmowania dwóch bardzo podobnych liczb. Korzystniejsza numerycznie jest postać uzyskana przez usunięcie niewymierności z licznika:

\[\sqrt{x+1}-\sqrt{x}=\frac{1}{\sqrt{x+1}+\sqrt{x}}.\]

Oba wzory są matematycznie równoważne, lecz drugi jest zwykle stabilniejszy numerycznie.

Porównywanie liczb

Wynik obliczeń zmiennoprzecinkowych nie zawsze powinien być porównywany z wartością oczekiwaną za pomocą dokładnej równości.

Zamiast sprawdzać:

\[x=y,\]

często sprawdza się, czy różnica jest dostatecznie mała:

\[|x-y|\leq \varepsilon.\]

Dla liczb o bardzo różnym rzędzie wielkości lepiej stosować jednocześnie tolerancję bezwzględną i względną, na przykład:

\[|x-y|\leq \varepsilon_{\mathrm{abs}}+\varepsilon_{\mathrm{rel}}\max(|x|,|y|).\]

Wartość tolerancji należy dobrać do skali danych, dokładności pomiarów i rodzaju wykonywanych obliczeń. Nie powinna być ona automatycznie utożsamiana z epsilonem maszynowym.

Sumowanie wielu składników

Podczas sumowania bardzo dużej liczby wartości niewielkie błędy mogą się kumulować. Znaczenie ma również kolejność składników.

Dokładność można poprawić między innymi przez:

• sumowanie liczb od najmniejszego do największego modułu,
• sumowanie parami w strukturze drzewiastej,
• zastosowanie algorytmu sumowania Kahana,
• wykorzystanie większej precyzji obliczeń.

Błąd bezwzględny i względny

Jeżeli dokładna wartość wynosi \(x\), a wartość zapisana lub obliczona przez komputer jest równa \(\operatorname{fl}(x)\), błąd bezwzględny wynosi

\[\Delta x=|x-\operatorname{fl}(x)|.\]

Błąd względny określa stosunek błędu bezwzględnego do modułu wartości dokładnej:

\[\delta x=\frac{|x-\operatorname{fl}(x)|}{|x|}.\]

Błąd względny pozwala ocenić dokładność niezależnie od skali liczby. Nie można go jednak bezpośrednio stosować dla \(x=0\), ponieważ prowadziłoby to do dzielenia przez zero.

Arytmetyka zmiennoprzecinkowa a metody numeryczne

W metodach numerycznych wynik zależy nie tylko od błędu wynikającego z przybliżonej metody, lecz również od błędów powstających podczas wykonywania działań zmiennoprzecinkowych.

Można wyróżnić między innymi:

• błąd danych — wynikający z niedokładności wartości wejściowych,
• błąd metody — wynikający z zastąpienia dokładnego problemu jego przybliżeniem,
• błąd obcięcia — powstający po przerwaniu nieskończonego procesu,
• błąd zaokrąglenia — wynikający z ograniczonej precyzji reprezentacji liczb.

Zmniejszenie kroku obliczeniowego nie zawsze poprawia wynik bez ograniczeń. Mniejszy krok może zmniejszyć błąd metody, ale jednocześnie zwiększyć liczbę działań i udział błędów zaokrągleń.

W praktyce często istnieje optymalny zakres wartości kroku, w którym łączny błąd jest najmniejszy.

Zjawiska te mają znaczenie między innymi przy różniczkowaniu numerycznym, całkowaniu, rozwiązywaniu równań nieliniowych oraz obliczeniach iteracyjnych. Dotyczą także dopasowywania funkcji i przetwarzania danych, takich jak zagadnienia omówione w artykule Aproksymacja i interpolacja.

Jak ograniczać błędy obliczeń?

• Stosować precyzję odpowiednią do problemu.
• Unikać odejmowania prawie równych liczb, gdy możliwe jest przekształcenie wzoru.
• Nie porównywać bezpośrednio wyników przybliżonych za pomocą dokładnej równości.
• Kontrolować rząd wielkości wartości pośrednich.
• Unikać niepotrzebnego mnożenia bardzo dużych i bardzo małych liczb.
• Skalować dane przed rozpoczęciem obliczeń.
• Stosować algorytmy numerycznie stabilne.
• Analizować wpływ zmiany kroku, tolerancji i precyzji na wynik.
• Sprawdzać występowanie wartości NaN i nieskończoności.
• Porównywać wynik z rozwiązaniem wzorcowym lub obliczeniami wykonanymi większą precyzją.

Większa liczba wykonanych działań nie musi oznaczać większej dokładności. Każde działanie zmiennoprzecinkowe może wprowadzać kolejny błąd zaokrąglenia.

Podsumowanie

Arytmetyka zmiennoprzecinkowa umożliwia zapisywanie liczb o bardzo różnym rzędzie wielkości za pomocą ograniczonej liczby bitów. Liczba jest reprezentowana przez znak, cechę i część znaczącą, tradycyjnie nazywaną mantysą.

Standard IEEE 754 definiuje formaty pojedynczej i podwójnej precyzji, zasady zaokrąglania oraz wartości specjalne, takie jak zero ze znakiem, liczby subnormalne, nieskończoność i NaN.

Nie wszystkie ułamki dziesiętne można dokładnie przedstawić w systemie binarnym. Liczby takie jak \(0{,}1\) i \(0{,}2\) są zapisywane jako skończone przybliżenia, co może prowadzić do niewielkich różnic w wynikach obliczeń.

Ważnym parametrem formatu jest epsilon maszynowy, określający odstęp między jedynką a następną większą liczbą reprezentowalną. Nie należy go jednak mylić z najmniejszą dodatnią liczbą ani automatycznie stosować jako tolerancji we wszystkich porównaniach.

Świadomość ograniczeń arytmetyki zmiennoprzecinkowej pozwala dobierać stabilniejsze algorytmy, właściwie interpretować wyniki i unikać wielu trudnych do wykrycia błędów obliczeniowych.

Utworzono: 12.06.2026 | Zmodyfikowano: 25.07.2026

Metody numeryczne – czym są i do czego służą?

Mgr inż.

Sebastian Dziarmaga-Działyński

Wirtualna kancelaria korepetytorska i konsultacyjna – wszechwiedza.pl

Metody numeryczne pozwalają rozwiązywać problemy matematyczne za pomocą skończonej liczby działań arytmetycznych wykonywanych ręcznie lub przez komputer. Stosuje się je przede wszystkim wtedy, gdy dokładne rozwiązanie analityczne jest trudne do otrzymania, nie istnieje w postaci elementarnej albo jego wyznaczenie wymagałoby zbyt wielu obliczeń.

Wynik metody numerycznej jest najczęściej wynikiem przybliżonym. Nie oznacza to jednak, że jest on przypadkowy lub pozbawiony wartości. Odpowiednio dobrany algorytm umożliwia kontrolowanie błędu, ocenę dokładności rozwiązania oraz osiągnięcie wyniku wystarczającego do zastosowań naukowych, technicznych i praktycznych.

Czym zajmują się metody numeryczne?

Metody numeryczne są działem matematyki stosowanej zajmującym się opracowywaniem i analizą sposobów otrzymywania liczbowych rozwiązań problemów matematycznych. Zamiast poszukiwać wyłącznie dokładnego wzoru opisującego rozwiązanie, tworzy się algorytm prowadzący do wartości przybliżonej.

Algorytmy numeryczne wykorzystują między innymi podstawowe działania arytmetyczne, porównania, powtarzanie określonych kroków oraz sprawdzanie warunków zakończenia obliczeń. Dzięki temu mogą być realizowane przez komputery, kalkulatory, sterowniki przemysłowe oraz inne urządzenia cyfrowe.

Metody numeryczne są używane na przykład do:

• wyznaczania miejsc zerowych funkcji i rozwiązywania równań nieliniowych,
• rozwiązywania układów równań liniowych,
• przybliżania funkcji na podstawie danych pomiarowych,
• obliczania wartości całek oznaczonych,
• przybliżonego wyznaczania pochodnych,
• rozwiązywania równań różniczkowych,
• wyznaczania wartości i wektorów własnych macierzy,
• optymalizacji funkcji,
• prowadzenia symulacji komputerowych.

Rozwiązanie analityczne i rozwiązanie numeryczne

W matematyce często rozróżnia się rozwiązania analityczne i numeryczne. Rozwiązanie analityczne ma zazwyczaj postać dokładnego wzoru. Rozwiązanie numeryczne jest natomiast określoną wartością liczbową lub zbiorem wartości, które z odpowiednią dokładnością przybliżają poszukiwane rozwiązanie.

Proste równanie liniowe można bez trudu rozwiązać analitycznie. Na przykład z równania

\[2x+4=10\]

otrzymujemy dokładnie \(x=3\). Nie każde równanie daje się jednak rozwiązać w tak prosty sposób. Przykładem może być równanie

\[x=\cos x.\]

Jego rozwiązania nie można przedstawić za pomocą prostego wzoru złożonego ze standardowych funkcji elementarnych. Można jednak zastosować odpowiednią metodę iteracyjną i otrzymać wartość przybliżoną

\[x\approx 0{,}739085.\]

Podobna sytuacja występuje przy obliczaniu niektórych całek oznaczonych. Jeżeli funkcji pierwotnej nie da się zapisać za pomocą funkcji elementarnych albo dostępne są jedynie dane tabelaryczne, wartość całki można wyznaczyć za pomocą metody prostokątów, trapezów, Simpsona lub bardziej zaawansowanych kwadratur numerycznych.

Algorytm numeryczny

Podstawą każdej metody numerycznej jest algorytm, czyli jednoznaczny przepis określający kolejność wykonywania działań. Algorytm otrzymuje dane wejściowe, przetwarza je zgodnie z ustalonymi regułami, a następnie zwraca wynik.

Typowy algorytm numeryczny może składać się z następujących etapów:

1. wprowadzenie danych i określenie poszukiwanego rozwiązania,
2. wybór metody numerycznej,
3. przyjęcie wartości początkowych lub początkowego przybliżenia,
4. wykonanie kolejnych działań lub iteracji,
5. oszacowanie błędu otrzymanego wyniku,
6. sprawdzenie warunku zakończenia obliczeń,
7. zwrócenie wyniku z określoną dokładnością.

Niektóre algorytmy prowadzą do rozwiązania po z góry ustalonej liczbie działań. Inne mają charakter iteracyjny. W metodzie iteracyjnej wyznacza się ciąg kolejnych przybliżeń

\[x_0,\ x_1,\ x_2,\ldots,\ x_n,\]

które — przy odpowiednich założeniach — zbliżają się do poszukiwanego rozwiązania. Obliczenia są przerywane, gdy różnica między kolejnymi przybliżeniami, wartość błędu lub inna miara dokładności spadnie poniżej ustalonej tolerancji.

Metoda numeryczna nie polega jedynie na otrzymaniu liczby. Równie ważne jest określenie, z jaką dokładnością liczba ta przybliża rzeczywiste rozwiązanie.

Dlaczego potrzebne są rozwiązania przybliżone?

W wielu rzeczywistych problemach rozwiązanie dokładne jest nieosiągalne lub niepraktyczne. Przyczyną może być złożoność modelu matematycznego, duża liczba niewiadomych, brak wzoru analitycznego albo niedokładność danych wejściowych.

Przykładowo analiza rozkładu temperatury w budynku, naprężeń w konstrukcji, przepływu powietrza wokół pojazdu lub zmian pogody wymaga rozwiązania bardzo rozbudowanych układów równań. Nawet jeżeli równania te wynikają z dobrze znanych praw fizyki, ich dokładne rozwiązanie może być niemożliwe.

W takich przypadkach badany obszar dzieli się na skończoną liczbę mniejszych elementów, a proces ciągły zastępuje się modelem dyskretnym. Komputer wykonuje następnie dużą liczbę działań i wyznacza przybliżone wartości poszukiwanych wielkości w wybranych punktach.

Dokładność, zbieżność i stabilność

Do najważniejszych pojęć używanych przy ocenie metod numerycznych należą dokładność, zbieżność, stabilność oraz koszt obliczeniowy.

Dokładność

Dokładność określa, jak blisko wartości dokładnej znajduje się wynik numeryczny. Jeżeli wartość dokładna wynosi \(x\), a wartość przybliżona \(\widetilde{x}\), można analizować między innymi błąd bezwzględny

\[\Delta x=\left|x-\widetilde{x}\right|\]

oraz błąd względny

\[\delta x=\frac{\left|x-\widetilde{x}\right|}{|x|}.\]

W praktyce wartość dokładna często nie jest znana. Stosuje się wówczas oszacowania błędu wynikające z właściwości metody albo porównuje się kolejne przybliżenia.

Zbieżność

Metoda jest zbieżna, jeżeli wraz ze zwiększaniem liczby iteracji, zmniejszaniem kroku obliczeń lub zagęszczaniem podziału jej wynik zbliża się do rozwiązania dokładnego.

Sama możliwość wykonywania kolejnych iteracji nie gwarantuje zbieżności. Źle dobrana metoda albo niewłaściwe przybliżenie początkowe mogą sprawić, że kolejne wartości będą oddalały się od rozwiązania lub zaczną oscylować.

Stabilność numeryczna

Stabilność opisuje zachowanie algorytmu w obecności niewielkich błędów danych i zaokrągleń. Stabilny algorytm nie powinien nadmiernie wzmacniać drobnych zaburzeń pojawiających się w trakcie obliczeń.

Należy przy tym odróżnić stabilność algorytmu od uwarunkowania samego problemu. Zadanie źle uwarunkowane może reagować dużą zmianą wyniku nawet na bardzo niewielką zmianę danych, niezależnie od jakości zastosowanego algorytmu.

Koszt obliczeniowy

Dwie metody mogą prowadzić do wyników o podobnej dokładności, lecz wymagać różnej liczby działań, ilości pamięci lub czasu pracy komputera. W praktyce poszukuje się więc nie tylko metod dokładnych i stabilnych, ale również wystarczająco wydajnych.

Skąd biorą się błędy obliczeń numerycznych?

Błędy numeryczne mogą pojawiać się na różnych etapach rozwiązywania problemu. Nie zawsze są skutkiem pomyłki człowieka lub wadliwego programu. Część z nich wynika z samej natury obliczeń przybliżonych.

• Błędy danych wejściowych wynikają na przykład z ograniczonej dokładności pomiaru.
• Błędy modelu powstają wtedy, gdy rzeczywiste zjawisko zostaje opisane za pomocą uproszczonych równań.
• Błędy metody są skutkiem zastąpienia dokładnego problemu jego przybliżeniem.
• Błędy obcięcia pojawiają się po przerwaniu nieskończonego procesu po skończonej liczbie kroków.
• Błędy zaokrągleń wynikają ze skończonej liczby cyfr używanych do zapisywania liczb w komputerze.
• Błędy implementacji mogą być spowodowane niepoprawnym programem lub niewłaściwym zastosowaniem algorytmu.

Przykładem ograniczeń arytmetyki komputerowej jest pozornie proste działanie

\[0{,}1+0{,}2.\]

W wielu językach programowania wynik zapisany z dużą liczbą cyfr może mieć postać zbliżoną do

\[0{,}30000000000000004.\]

Nie oznacza to błędu działania procesora. Liczby \(0{,}1\) i \(0{,}2\) nie mają skończonego rozwinięcia w systemie dwójkowym, dlatego muszą zostać zapisane w pamięci komputera w sposób przybliżony. Zagadnienia te są związane z arytmetyką zmiennoprzecinkową i pojęciem epsilonu maszynowego.

Główne dziedziny metod numerycznych

Metody numeryczne obejmują wiele grup algorytmów. Poszczególne grupy odpowiadają różnym typom problemów matematycznych.

Rozwiązywanie równań nieliniowych

Celem jest wyznaczenie takich wartości \(x\), dla których

\[f(x)=0.\]

Do podstawowych metod należą metoda bisekcji, metoda regula falsi, metoda siecznych i metoda Newtona. Różnią się one szybkością zbieżności, wymaganiami dotyczącymi funkcji oraz sposobem wyboru wartości początkowych.

Układy równań liniowych

Wiele problemów technicznych prowadzi do układu zapisanego w postaci macierzowej

\[A\mathbf{x}=\mathbf{b}.\]

Układy takie rozwiązuje się między innymi za pomocą eliminacji Gaussa, rozkładu LU oraz metod iteracyjnych Jacobiego i Gaussa – Seidla. Podstawowe informacje o działaniach wykonywanych na macierzach przedstawiono w artykule Macierze – część 2.

Interpolacja i aproksymacja

Interpolacja polega na wyznaczeniu funkcji przechodzącej przez zadane punkty. Aproksymacja ma natomiast na celu znalezienie funkcji, która możliwie dobrze przybliża dane, ale nie musi przechodzić dokładnie przez każdy punkt.

Zagadnienia interpolacji liniowej, wielomianów Lagrange’a i Newtona oraz aproksymacji średniokwadratowej zostały szerzej omówione w artykule Aproksymacja i interpolacja.

Różniczkowanie numeryczne

Jeżeli funkcja jest znana jedynie w wybranych punktach, jej pochodną można przybliżyć za pomocą ilorazów różnicowych. Przykładowo pochodną w punkcie \(x\) można oszacować wzorem

\[f'(x)\approx\frac{f(x+h)-f(x)}{h}.\]

Dokładniejsze przybliżenia można otrzymać za pomocą różnic centralnych. Wprowadzenie do pojęcia pochodnej znajduje się w artykule Pochodna funkcji.

Całkowanie numeryczne

Całkowanie numeryczne służy do przybliżonego wyznaczania wartości całki oznaczonej

\[I=\int_a^b f(x)\,dx.\]

Przedział całkowania dzieli się na mniejsze części, a pole pod wykresem funkcji zastępuje się sumą pól prostszych figur. W metodzie prostokątów są to prostokąty, w metodzie trapezów — trapezy, natomiast metoda Simpsona wykorzystuje przybliżenie funkcji wielomianami drugiego stopnia.

Całkowanie numeryczne jest szczególnie przydatne wtedy, gdy nie można znaleźć funkcji pierwotnej w postaci elementarnej albo gdy wartości funkcji pochodzą z pomiarów. Analityczne metody obliczania całek zostały przedstawione między innymi w artykule Całki nieoznaczone – metody podstawowe.

Równania różniczkowe

Równania różniczkowe opisują zależności między funkcją i jej pochodnymi. Są wykorzystywane do modelowania ruchu, przepływu ciepła, drgań, obwodów elektrycznych, wzrostu populacji oraz wielu innych procesów.

Jeżeli równanie ma postać

\[y'=f(x,y),\qquad y(x_0)=y_0,\]

najprostszą metodą numeryczną jest metoda Eulera:

\[y_{n+1}=y_n+h f(x_n,y_n).\]

Dokładniejsze wyniki można otrzymać za pomocą ulepszonej metody Eulera oraz metod Rungego – Kutty. Ważnym zagadnieniem jest przy tym odpowiedni dobór kroku \(h\).

Zastosowania metod numerycznych

Metody numeryczne są obecne niemal wszędzie tam, gdzie matematyczny model rzeczywistości musi zostać przekształcony w konkretny wynik liczbowy.

• Mechanika i budownictwo: analiza naprężeń, odkształceń, drgań i stateczności konstrukcji.
• Elektrotechnika: obliczanie prądów i napięć w rozbudowanych obwodach.
• Energetyka: modelowanie przepływu ciepła, pracy sieci energetycznych i procesów spalania.
• Mechanika płynów: symulowanie przepływu cieczy i gazów.
• Astronomia: wyznaczanie orbit i przewidywanie położeń ciał niebieskich.
• Ekonomia i finanse: optymalizacja, wycena instrumentów finansowych i prognozowanie.
• Geodezja: wyrównywanie obserwacji i przetwarzanie wyników pomiarów.
• Grafika komputerowa: modelowanie krzywych, powierzchni, oświetlenia i ruchu.
• Uczenie maszynowe: numeryczna optymalizacja parametrów modeli.
• Medycyna: rekonstrukcja obrazów, modelowanie przepływu krwi i analiza sygnałów.

Czy komputer zawsze podaje poprawny wynik?

Komputer wykonuje niezwykle szybko ogromną liczbę działań, ale nie ocenia samodzielnie matematycznego sensu otrzymanego wyniku. Poprawnie napisany program może dokładnie zrealizować algorytm, który został niewłaściwie dobrany do rozwiązywanego problemu.

Wynik może być nieprawidłowy lub niewiarygodny między innymi wtedy, gdy:

• dane wejściowe zawierają duże błędy,
• problem jest źle uwarunkowany,
• metoda nie jest zbieżna dla wybranych danych,
• przyjęto zbyt duży krok obliczeniowy,
• iteracje przerwano zbyt wcześnie,
• nastąpiła utrata cyfr znaczących,
• wynik przekroczył zakres reprezentowanych liczb,
• program zawiera błąd implementacyjny.

Dlatego wynik obliczeń numerycznych powinien być zawsze analizowany. Warto sprawdzić jego jednostkę, rząd wielkości, zgodność z warunkami zadania oraz zachowanie po niewielkiej zmianie parametrów obliczeń.

Metody numeryczne a rozwój komputerów

Metody przybliżone były stosowane jeszcze przed powstaniem komputerów elektronicznych. Obliczenia wykonywano ręcznie, za pomocą tablic matematycznych, suwaków logarytmicznych i mechanicznych maszyn liczących. Ze względu na czasochłonność możliwe było jednak rozwiązywanie jedynie stosunkowo niewielkich problemów.

Rozwój komputerów umożliwił realizowanie algorytmów wymagających milionów, miliardów, a obecnie nawet znacznie większej liczby operacji. Dzięki temu metody numeryczne stały się jednym z podstawowych narzędzi współczesnej nauki, techniki i analizy danych.

Wzrost mocy obliczeniowej nie eliminuje jednak potrzeby badania dokładności i stabilności. Szybszy komputer wykona niestabilny lub błędny algorytm szybciej, ale nie sprawi automatycznie, że otrzymany wynik będzie poprawny.

Podsumowanie

Metody numeryczne umożliwiają otrzymywanie przybliżonych rozwiązań problemów, których nie można wygodnie rozwiązać analitycznie. Ich podstawą są algorytmy realizowane za pomocą skończonej liczby działań i prowadzące do wyniku o kontrolowanej dokładności.

Przy wyborze metody należy uwzględnić nie tylko dokładność, ale również zbieżność, stabilność i koszt obliczeniowy. Ważna jest także analiza źródeł błędów, w tym błędów danych, modelu, obcięcia i zaokrągleń.

Do najważniejszych dziedzin metod numerycznych należą rozwiązywanie równań i układów równań, interpolacja i aproksymacja, różniczkowanie i całkowanie numeryczne, rozwiązywanie równań różniczkowych oraz optymalizacja. Zagadnienia te tworzą podstawę obliczeń komputerowych stosowanych we współczesnej nauce i technice.

Utworzono: 10.06.2026 | Zmodyfikowano: 24.07.2026

Krzywa Lorenza i współczynnik Giniego — miary koncentracji

Mgr inż.

Sebastian Dziarmaga-Działyński

Wirtualna kancelaria korepetytorska i konsultacyjna – wszechwiedza.pl

Krzywa Lorenza i współczynnik Giniego służą do badania koncentracji wartości cechy. Pokazują, czy suma badanej cechy jest rozłożona równomiernie między jednostki, czy też skupia się w rękach niewielkiej części zbiorowości.

W statystyce opisowej mówimy o różnych rodzajach miar. Miary położenia pokazują, gdzie znajduje się typowy poziom danych. Miary zróżnicowania opisują rozproszenie wartości. Miary asymetrii informują, czy rozkład ma dłuższy ogon po lewej lub prawej stronie. Z kolei kurtoza i eksces opisują koncentrację rozkładu wokół środka oraz zachowanie ogonów.

Krzywa Lorenza i współczynnik Giniego należą do innego ujęcia koncentracji. Nie pytają o to, czy rozkład jest spiczasty albo spłaszczony, lecz o to, jak nierównomiernie suma wartości cechy jest rozdzielona między jednostki. Najczęściej stosuje się je do analizy dochodów, majątku, sprzedaży, udziałów rynkowych albo innych wielkości nieujemnych, które można sumować.

Uwaga językowa: poprawna nazwa to krzywa Lorenza, a nie „krzywa Lorentza”. Współczynnik zapisujemy jako współczynnik Giniego.

Czym jest koncentracja wartości cechy?

Koncentracja wartości cechy oznacza, że duża część łącznej sumy tej cechy przypada na niewielką część jednostek. Przykładowo, jeżeli w pewnej grupie kilka osób posiada większość całego majątku, mówimy o silnej koncentracji majątku.

Podobnie może być ze sprzedażą przedsiębiorstw, liczbą klientów, powierzchnią gospodarstw rolnych albo udziałami w rynku. W każdym z tych przypadków interesuje nas nie tylko to, jakie wartości mają poszczególne jednostki, ale także to, jaka część całkowitej sumy przypada na kolejne części zbiorowości.

Współczynnik Giniego i krzywa Lorenza są szczególnie użyteczne wtedy, gdy badana cecha jest nieujemna i addytywna, czyli można sensownie sumować jej wartości. Dlatego dobrze nadają się do analizy dochodów, majątku, produkcji, sprzedaży albo udziałów rynkowych.

Krzywa Lorenza — intuicja

Krzywa Lorenza pokazuje, jaka część łącznej sumy cechy przypada na kolejne części zbiorowości uporządkowanej od najmniejszych wartości do największych.

Na osi poziomej odkładamy skumulowany udział jednostek, na przykład osób, gospodarstw domowych albo przedsiębiorstw. Na osi pionowej odkładamy skumulowany udział sumy cechy, na przykład dochodów, majątku albo sprzedaży.

Jeżeli analizujemy dochody, to krzywa Lorenza odpowiada na pytania w rodzaju:

• jaki procent łącznego dochodu przypada na najbiedniejsze \(25\%\) osób,
• jaki procent łącznego dochodu przypada na najbiedniejsze \(50\%\) osób,
• jaki procent łącznego dochodu przypada na najbiedniejsze \(75\%\) osób,
• jak bardzo rzeczywisty rozkład dochodów odbiega od pełnej równości.

Przekątna pełnej równości

Na wykresie krzywej Lorenza bardzo ważna jest przekątna pełnej równości. Jest to linia, która biegnie od punktu \((0,0)\) do punktu \((1,1)\), czyli od \(0\%\) jednostek i \(0\%\) sumy cechy do \(100\%\) jednostek i \(100\%\) sumy cechy.

Gdyby każdy miał dokładnie tyle samo, krzywa Lorenza pokrywałaby się z tą przekątną. Wtedy:

• \(25\%\) osób miałoby \(25\%\) łącznego dochodu,
• \(50\%\) osób miałoby \(50\%\) łącznego dochodu,
• \(75\%\) osób miałoby \(75\%\) łącznego dochodu,
• \(100\%\) osób miałoby \(100\%\) łącznego dochodu.

W rzeczywistych danych krzywa Lorenza zwykle przebiega poniżej przekątnej pełnej równości. Im bardziej oddala się od przekątnej, tym większa jest koncentracja wartości cechy.

Jak zbudować krzywą Lorenza?

Aby zbudować krzywą Lorenza, wykonujemy kilka kroków:

1. Porządkujemy jednostki rosnąco według wartości badanej cechy.
2. Obliczamy sumę wszystkich wartości cechy.
3. Obliczamy skumulowany udział jednostek.
4. Obliczamy skumulowany udział wartości cechy.
5. Na wykresie zaznaczamy punkty i łączymy je linią.

Dla uporządkowanych wartości:

\[ x_1 \leq x_2 \leq \ldots \leq x_n \]

skumulowany udział jednostek dla \(i\)-tej pozycji wynosi:

\[ P_i=\frac{i}{n} \]

a skumulowany udział wartości cechy:

\[ Q_i=\frac{x_1+x_2+\ldots+x_i}{x_1+x_2+\ldots+x_n} \]

Dodatkowo przyjmujemy punkt początkowy:

\[ P_0=0,\qquad Q_0=0 \]

Prosty przykład krzywej Lorenza

Załóżmy, że cztery osoby osiągają dochody:

\[ 2,\ 2,\ 6,\ 10 \]

Dane są już uporządkowane rosnąco. Suma dochodów wynosi:

\[ 2+2+6+10=20 \]

Budujemy tabelę udziałów skumulowanych:

Interpretacja jest następująca:

• najbiedniejsze \(25\%\) osób ma \(10\%\) łącznego dochodu,
• najbiedniejsze \(50\%\) osób ma \(20\%\) łącznego dochodu,
• najbiedniejsze \(75\%\) osób ma \(50\%\) łącznego dochodu,
• wszystkie osoby razem mają \(100\%\) dochodu.

Gdyby dochody były rozłożone idealnie równo, to \(50\%\) osób miałoby \(50\%\) dochodu. W naszym przykładzie \(50\%\) osób ma tylko \(20\%\) dochodu, więc widzimy wyraźną nierównomierność.

Współczynnik Giniego

Współczynnik Giniego jest liczbową miarą koncentracji wartości cechy. Pokazuje, jak daleko rzeczywisty rozkład odbiega od pełnej równości.

Dla nieujemnych wartości cechy współczynnik Giniego przyjmuje zwykle wartości od \(0\) do \(1\). Czasem podaje się go także w procentach, od \(0\%\) do \(100\%\).

Warto pamiętać, że współczynnik Giniego nie mówi, czy poziom dochodów jest wysoki czy niski. Dwie grupy mogą mieć taki sam współczynnik Giniego, ale zupełnie inny przeciętny poziom dochodów. Dlatego współczynnik Giniego warto interpretować razem z innymi miarami, na przykład średnią, medianą i kwartylami.

Geometryczna interpretacja współczynnika Giniego

Współczynnik Giniego można interpretować geometrycznie na wykresie krzywej Lorenza. Niech:

• \(A\) oznacza pole między przekątną pełnej równości a krzywą Lorenza,
• \(B\) oznacza pole pod krzywą Lorenza.

Wtedy:

\[ G=\frac{A}{A+B} \]

Ponieważ całe pole pod przekątną pełnej równości wynosi \(\frac{1}{2}\), można też zapisać:

\[ G=2A \]

albo:

\[ G=1-2B \]

Ta interpretacja dobrze pokazuje sens miary: im bardziej krzywa Lorenza oddala się od przekątnej pełnej równości, tym większe pole \(A\), a więc tym większy współczynnik Giniego.

Wzór na współczynnik Giniego dla danych szczegółowych

Jednym z wygodnych wzorów na współczynnik Giniego dla uporządkowanych rosnąco, nieujemnych wartości:

\[ x_1 \leq x_2 \leq \ldots \leq x_n \]

jest wzór:

gdzie:

• \(x_i\) — uporządkowane rosnąco wartości cechy,
• \(i\) — numer pozycji w uporządkowanym szeregu,
• \(n\) — liczba jednostek.

Ten wzór jest wygodny przy niewielkich zbiorach danych, ponieważ pozwala obliczyć współczynnik Giniego bez rysowania krzywej Lorenza.

Przykład obliczania współczynnika Giniego

Obliczmy współczynnik Giniego dla danych:

\[ 2,\ 2,\ 6,\ 10 \]

Dane są już uporządkowane rosnąco. Najpierw obliczamy sumę wartości:

\[ \sum_{i=1}^{n}x_i = 2+2+6+10=20 \]

Następnie obliczamy sumę \(i x_i\):

\[ \sum_{i=1}^{n}i x_i =1\cdot 2+2\cdot 2+3\cdot 6+4\cdot 10 \]

\[ =2+4+18+40=64 \]

Podstawiamy do wzoru:

\[ G=\frac{2\cdot 64}{4\cdot 20}-\frac{4+1}{4} \]

\[ G=\frac{128}{80}-\frac{5}{4} \]

\[ G=1{,}6-1{,}25=0{,}35 \]

Współczynnik Giniego wynosi:

\[ G=0{,}35 \]

Oznacza to umiarkowaną koncentrację wartości cechy. Rozkład nie jest równy, ale nie mamy też do czynienia z sytuacją skrajnej koncentracji.

Przykład pełnej równości i silnej koncentracji

Dla porównania rozważmy trzy zbiory danych:

W zbiorze A każda jednostka ma taką samą wartość, więc współczynnik Giniego wynosi \(0\). W zbiorze B wartości są nierówne, ale nie cała suma przypada jednej jednostce. W zbiorze C prawie cała suma wartości cechy skupia się w jednej jednostce, dlatego koncentracja jest bardzo silna.

Współczynnik Giniego a średnia, mediana i asymetria

Współczynnik Giniego nie zastępuje średniej, mediany ani miar asymetrii. Odpowiada na inne pytanie.

Przykładowo dwa kraje mogą mieć podobny średni dochód, ale różny współczynnik Giniego. Oznaczałoby to, że przeciętny poziom dochodu jest podobny, ale nierówności dochodowe są inne.

Podobnie rozkład dochodów może być prawostronnie asymetryczny, ponieważ ma długi ogon wysokich wartości, ale współczynnik Giniego dodatkowo pokazuje, jak silnie łączna suma dochodu koncentruje się wśród bogatszej części populacji.

Kiedy warto stosować współczynnik Giniego?

Współczynnik Giniego warto stosować wtedy, gdy interesuje nas nierównomierność podziału sumy wartości cechy. Typowe zastosowania to:

• analiza nierówności dochodowych,
• analiza nierówności majątkowych,
• badanie koncentracji sprzedaży między przedsiębiorstwami,
• badanie koncentracji udziałów rynkowych,
• analiza koncentracji produkcji, powierzchni gospodarstw albo liczby klientów.

Nie należy jednak stosować współczynnika Giniego mechanicznie do każdej zmiennej. Najlepiej sprawdza się dla wartości nieujemnych, które można sumować i interpretować jako udziały w całości.

Ograniczenia współczynnika Giniego

Współczynnik Giniego jest bardzo użyteczny, ale ma też ograniczenia.

• Nie pokazuje, gdzie dokładnie występują nierówności — dwie różne krzywe Lorenza mogą dać podobną wartość współczynnika.
• Nie informuje o przeciętnym poziomie cechy — mówi o nierównomierności, a nie o wysokości dochodów czy majątku.
• Wymaga ostrożności przy danych z wartościami ujemnymi, ponieważ klasyczna interpretacja udziałów w sumie może wtedy tracić sens.
• Nie mówi nic o przyczynach nierówności — jest miarą opisową, a nie wyjaśniającą.

Najczęstsze błędy przy interpretacji współczynnika Giniego

1. Mylenie nierówności z poziomem dochodu

Współczynnik Giniego nie mówi, czy ludzie są bogaci czy biedni. Mówi tylko, jak nierównomiernie rozłożona jest suma dochodu. Możliwa jest sytuacja, w której dwa kraje mają podobny współczynnik Giniego, ale zupełnie inny poziom przeciętnych dochodów.

2. Pomijanie krzywej Lorenza

Sam współczynnik Giniego jest jedną liczbą. Krzywa Lorenza pokazuje więcej szczegółów, ponieważ pozwala zobaczyć, w której części rozkładu pojawia się największe odchylenie od pełnej równości.

3. Stosowanie miary do danych, dla których suma nie ma sensu

Współczynnik Giniego najlepiej działa dla cech, których wartości są nieujemne i mogą być sumowane. Nie ma sensu stosować go do zmiennych jakościowych ani do wielu skal, w których suma wartości nie ma naturalnej interpretacji.

4. Zapominanie o uporządkowaniu danych

Przy korzystaniu ze wzoru dla danych szczegółowych wartości należy uporządkować rosnąco. To samo dotyczy budowy krzywej Lorenza — zaczynamy od jednostek o najmniejszych wartościach cechy i przechodzimy do coraz większych.

Zadanie dla czytelnika

Podsumowanie

Krzywa Lorenza i współczynnik Giniego służą do badania koncentracji wartości cechy. Krzywa Lorenza pokazuje graficznie, jaka część sumy cechy przypada na kolejne części uporządkowanej zbiorowości. Przekątna pełnej równości pokazuje sytuację idealnie równomiernego podziału.

Współczynnik Giniego jest liczbową miarą odchylenia od pełnej równości. Wartość \(0\) oznacza pełną równość, a wartości bliższe \(1\) oznaczają coraz silniejszą koncentrację. Miara ta jest szczególnie przydatna przy analizie dochodów, majątku, sprzedaży, udziałów rynkowych i innych nieujemnych wielkości, które można sumować.

Najważniejsze jest to, aby nie mylić współczynnika Giniego z miarami położenia, asymetrii czy kurtozy. Gini odpowiada na osobne pytanie: jak nierównomiernie łączna suma cechy jest podzielona między jednostki.

Utworzono: 05.06.2026 | Zmodyfikowano: 21.06.2026

Kurtoza, eksces i koncentracja w statystyce opisowej

Mgr inż.

Sebastian Dziarmaga-Działyński

Wirtualna kancelaria korepetytorska i konsultacyjna – wszechwiedza.pl

Kurtoza, eksces i miary koncentracji opisują kształt rozkładu danych: czy wartości silnie skupiają się wokół środka, czy rozkład jest bardziej spłaszczony, a także czy w ogonach rozkładu pojawiają się wartości skrajne.

W poprzednich artykułach omawialiśmy miary położenia, miary zróżnicowania oraz miary asymetrii. Każda z tych grup opisuje inną cechę rozkładu. Miary położenia mówią, gdzie znajduje się środek danych. Miary zróżnicowania informują, jak bardzo dane są rozproszone. Miary asymetrii pokazują, czy rozkład ma dłuższy ogon po jednej stronie.

Kurtoza i eksces pozwalają spojrzeć na rozkład jeszcze inaczej. Interesuje nas tutaj przede wszystkim koncentracja wartości wokół środka rozkładu oraz zachowanie ogonów rozkładu. W klasycznej dydaktyce mówi się często o rozkładach bardziej „spiczastych” i bardziej „spłaszczonych”, choć warto pamiętać, że kurtoza wiąże się także z ciężkością ogonów i wpływem wartości skrajnych.

Ten artykuł jest częścią szerszego cyklu o miarach statystycznych. Ogólne omówienie całego podziału znajduje się w artykule Miary statystyczne w statystyce opisowej.

Czym jest koncentracja rozkładu?

Koncentracja rozkładu opisuje, w jakim stopniu wartości skupiają się wokół środka rozkładu oraz jak zachowują się wartości oddalone od środka. Jeżeli wiele obserwacji znajduje się bardzo blisko wartości typowej, rozkład może sprawiać wrażenie bardziej skupionego. Jeżeli obserwacje są bardziej równomiernie rozłożone albo występują szerokie obszary wartości, rozkład może być bardziej spłaszczony.

W praktyce dydaktycznej koncentrację rozkładu często porównuje się z rozkładem normalnym. Rozkład normalny pełni wtedy rolę punktu odniesienia: można zapytać, czy badany rozkład jest bardziej spiczasty, podobny, czy bardziej spłaszczony od rozkładu normalnego.

Kurtoza a eksces — ważna uwaga terminologiczna

Przy omawianiu kurtozy trzeba uważać na terminologię, ponieważ w różnych książkach, programach i kursach można spotkać nieco inne nazewnictwo.

W klasycznym ujęciu kurtoza to iloraz czwartego momentu centralnego i czwartej potęgi odchylenia standardowego. Dla rozkładu normalnego taka kurtoza wynosi \(3\).

Eksces to kurtoza pomniejszona o \(3\). Dzięki temu dla rozkładu normalnego eksces wynosi \(0\). Taka postać jest wygodna, ponieważ wynik od razu pokazuje, czy rozkład jest bardziej leptokurtyczny, mezokurtyczny czy platykurtyczny.

Można jednak spotkać autorów i programy, które słowem „kurtoza” nazywają właśnie eksces, czyli wartość już pomniejszoną o \(3\). Dlatego przy interpretacji wyniku zawsze warto sprawdzić, jakiej definicji używa podręcznik, wykładowca albo program statystyczny.

Jeżeli program zwraca wartość dodatnią dla rozkładu bardziej spiczastego od normalnego, a wartość ujemną dla rozkładu bardziej spłaszczonego od normalnego, to najprawdopodobniej raportuje eksces, a nie klasyczną kurtozę równą \(3\) dla rozkładu normalnego.

Czwarty moment centralny

Podstawą klasycznej miary koncentracji jest czwarty moment centralny. Jest on podobny konstrukcyjnie do wariancji, ale zamiast drugich potęg odchyleń od średniej wykorzystuje czwarte potęgi.

Przypomnijmy: wariancja opiera się na kwadratach odchyleń od średniej. Czwarty moment centralny idzie dalej i podnosi te odchylenia do czwartej potęgi. Dlatego bardzo duże odchylenia od średniej mają szczególnie silny wpływ na jego wartość. Samą wariancję i odchylenie standardowe omawiamy szerzej w artykule Odchylenie standardowe i wariancja.

gdzie:

• \(m_4\) — czwarty moment centralny,
• \(x_i\) — kolejna wartość cechy,
• \(\bar{x}\) — średnia arytmetyczna,
• \(n\) — liczba obserwacji.

Czwarty moment centralny jest miarą bezwzględną. Jego wartość zależy od jednostki pomiaru i skali danych, dlatego samodzielnie bywa trudny do interpretacji. Z tego powodu zwykle przechodzi się do miary względnej, czyli kurtozy.

Kurtoza jako klasyczny współczynnik koncentracji

Klasyczny współczynnik koncentracji, nazywany kurtozą, otrzymujemy przez podzielenie czwartego momentu centralnego przez czwartą potęgę odchylenia standardowego. Równoważnie można powiedzieć, że dzielimy czwarty moment centralny przez kwadrat wariancji.

gdzie:

• \(\beta_2\) — kurtoza, czyli klasyczny współczynnik koncentracji,
• \(m_4\) — czwarty moment centralny,
• \(s\) — odchylenie standardowe,
• \(s^2\) — wariancja.

Dzięki podzieleniu przez \(s^4\) miara staje się względna, czyli nie zależy bezpośrednio od jednostki pomiaru. Możemy więc porównywać kształt różnych rozkładów, nawet jeśli dane są wyrażone w innych jednostkach.

Eksces

Eksces jest kurtozą pomniejszoną o \(3\). Odejmowanie trójki wynika z tego, że klasyczna kurtoza rozkładu normalnego wynosi \(3\). Po odjęciu \(3\) rozkład normalny ma wartość odniesienia równą \(0\).

Jeżeli:

• \(\gamma_2=0\), rozkład ma koncentrację podobną do normalnego,
• \(\gamma_2>0\), rozkład jest leptokurtyczny,
• \(\gamma_2<0\), rozkład jest platykurtyczny.

W praktyce eksces jest często wygodniejszy niż sama kurtoza, ponieważ wynik dodatni, ujemny albo równy zero od razu wskazuje, jak rozkład wypada względem rozkładu normalnego.

Rozkład leptokurtyczny, mezokurtyczny i platykurtyczny

Najczęściej w statystyce opisowej wyróżnia się trzy typy rozkładów ze względu na kurtozę lub eksces: leptokurtyczne, mezokurtyczne i platykurtyczne. Punktem odniesienia jest rozkład normalny, dla którego klasyczna kurtoza wynosi \(3\), a eksces \(0\).

Rozkład mezokurtyczny

Rozkład mezokurtyczny ma koncentrację podobną do rozkładu normalnego. W klasycznej konwencji jego kurtoza wynosi około \(3\), a eksces około \(0\).

Można go traktować jako punkt odniesienia. Nie oznacza to, że każdy rozkład mezokurtyczny musi być dokładnie normalny, lecz że jego koncentracja i ogony są podobne do rozkładu normalnego.

Rozkład leptokurtyczny

Rozkład leptokurtyczny ma dodatni eksces:

\[ \gamma_2>0 \]

W klasycznym opisie mówi się, że taki rozkład jest bardziej spiczasty niż normalny. Warto jednak pamiętać, że wysoka kurtoza oznacza również większy wpływ wartości oddalonych od średniej. Dlatego rozkłady leptokurtyczne często wiążą się z cięższymi ogonami i większą skłonnością do wartości skrajnych.

Przykładowo w finansach dodatni eksces może oznaczać, że bardzo duże dodatnie lub ujemne stopy zwrotu pojawiają się częściej, niż sugerowałby prosty model normalny.

Rozkład platykurtyczny

Rozkład platykurtyczny ma ujemny eksces:

\[ \gamma_2<0 \]

W klasycznej dydaktyce mówi się, że taki rozkład jest bardziej spłaszczony niż normalny. Oznacza to mniejszą koncentrację wokół środka albo lżejsze ogony w porównaniu z rozkładem normalnym.

Przykładem rozkładu o ujemnym ekscesie może być rozkład, w którym wartości są bardziej równomiernie rozłożone w pewnym przedziale, bez silnego skupienia w środku i bez bardzo długich ogonów.

Kurtoza nie jest tym samym co asymetria

Kurtoza i asymetria opisują dwie różne cechy rozkładu. Asymetria mówi, czy rozkład ma dłuższy ogon po lewej czy po prawej stronie. Kurtoza mówi o koncentracji wartości i zachowaniu ogonów rozkładu względem rozkładu normalnego.

Oznacza to, że rozkład może być jednocześnie:

• symetryczny i leptokurtyczny,
• symetryczny i platykurtyczny,
• prawostronnie asymetryczny i leptokurtyczny,
• prawostronnie asymetryczny i platykurtyczny,
• lewostronnie asymetryczny i leptokurtyczny,
• lewostronnie asymetryczny i platykurtyczny.

To bardzo ważne, ponieważ nie należy utożsamiać „spiczastości” rozkładu z jego asymetrią. Rozkład może być wyraźnie przechylony w prawo, a jednocześnie mieć małą koncentrację w środku. Może też być symetryczny, ale mieć bardzo ciężkie ogony. Miary asymetrii omawiamy osobno w artykule Miary asymetrii w statystyce opisowej.

Jak interpretować wynik kurtozy i ekscesu?

Przy interpretacji trzeba najpierw ustalić, czy mamy do czynienia z kurtozą, czy z ekscesem. To najczęstsze źródło nieporozumień.

Dlatego przed interpretacją warto sprawdzić, czy program lub podręcznik podaje kurtozę w postaci \(\beta_2\), czy eksces w postaci \(\gamma_2=\beta_2-3\).

Kurtoza w programach statystycznych i arkuszach kalkulacyjnych

W programach statystycznych i arkuszach kalkulacyjnych trzeba zachować szczególną ostrożność, ponieważ nazwa funkcji może sugerować „kurtozę”, ale wynik bywa raportowany względem rozkładu normalnego. W praktyce oznacza to, że funkcja może zwracać eksces, czyli wartość znormalizowaną tak, aby rozkład normalny miał wynik równy \(0\).

Jeżeli program zwraca wynik dodatni, oznacza to zwykle rozkład bardziej leptokurtyczny niż normalny. Jeżeli wynik jest ujemny, oznacza rozkład bardziej platykurtyczny niż normalny. Jeżeli wynik jest bliski zeru, rozkład jest pod względem kurtozy zbliżony do normalnego.

Najczęstsze błędy przy interpretacji kurtozy

1. Mylenie kurtozy z ekscesem

Najczęstszy błąd polega na tym, że ktoś otrzymuje wynik \(0\) i uznaje, że kurtoza wynosi zero. W wielu programach taki wynik oznacza jednak eksces równy zero, czyli kurtozę równą \(3\), a więc rozkład podobny do normalnego pod względem koncentracji.

2. Utożsamianie kurtozy wyłącznie z wysokością wierzchołka

W podstawowych wyjaśnieniach mówi się często, że rozkład leptokurtyczny jest bardziej spiczasty, a platykurtyczny bardziej spłaszczony. To użyteczna intuicja, ale niepełna. Kurtoza jest silnie związana także z zachowaniem ogonów i wartościami oddalonymi od średniej.

3. Mylenie koncentracji z asymetrią

Rozkład może być symetryczny, ale leptokurtyczny albo platykurtyczny. Może też być asymetryczny i jednocześnie mieć dodatni lub ujemny eksces. Asymetria i kurtoza opisują różne cechy rozkładu, dlatego należy je interpretować osobno.

4. Interpretowanie kurtozy bez sprawdzenia skali i wartości odstających

Pojedyncze wartości skrajne mogą mieć duży wpływ na czwarty moment centralny, a więc także na kurtozę. Wynika to z czwartej potęgi odchyleń od średniej. Dlatego przy wysokiej kurtozie warto sprawdzić, czy w danych nie występują obserwacje nietypowe.

Zadanie dla czytelnika

Podsumowanie

Kurtoza, eksces i miary koncentracji opisują kształt rozkładu pod względem skupienia wartości wokół środka oraz zachowania ogonów. Klasyczna kurtoza jest ilorazem czwartego momentu centralnego i czwartej potęgi odchylenia standardowego. Dla rozkładu normalnego wynosi \(3\).

Eksces to kurtoza pomniejszona o \(3\). Dzięki temu rozkład normalny ma eksces równy \(0\). Dodatni eksces oznacza rozkład leptokurtyczny, a ujemny eksces oznacza rozkład platykurtyczny.

Najważniejsze jest to, aby nie mylić kurtozy z ekscesem oraz nie utożsamiać koncentracji z asymetrią. Asymetria mówi o stronie ogona rozkładu, natomiast kurtoza i eksces opisują koncentrację oraz zachowanie ogonów względem rozkładu normalnego.

Utworzono: 05.06.2026 | Zmodyfikowano: 21.06.2026