microgeek.eu

Tematyka związana z Cyfrowym przetwarzaniem sygnałów siłą
rzeczy „dotyka” zagadnień numerycznych. To już samo w sobie
wykracza poza ramy samego DSP i każdy z tych elementów może
żyć własnym życiem. Jedynym elementem ich scalającym jest pewnego
rodzaju podobieństwo zastosowania. Z tego powodu zdecydowałem się
na utworzenie niezależnego wątku poświęconego szeroko pojętym
zagadnieniom numerycznym (rozwiązywania określonych problemów
matematycznych posiłkując się liczbami zmiennoprzecinkowymi
i językiem programowania → w tym konkretnym przypadku językiem C).
W sumie wszystko opiera się o wiedzę uniwersalną, wypracowaną przez
wysiłek wielu ludzi, chociaż w paru przypadkach kierowałem się własną
koncepcją oraz potrzebą. Z pewnością niektóre rozwiązania można
zrealizować inaczej (nie znaczy, że lepiej lub gorzej → po prostu inaczej).

Numeryczne poszukiwanie zer wielomianu

Potrzeba i koncepcja

W wielu zagadnieniach technicznych zachodzi potrzeba znalezienia pierwiastków wielomianów (rozwiązania równania z wielomianem). Zakładając, że wielomian stopnia n ma współczynniki rzeczywiste, to jego rozwiązanie będzie zawierać n rozwiązań (rzeczywistych oraz zespolonych). W przypadku rozwiązań zespolonych, będą one sprzężone.
Algorytm rozwiązania opiera się o metodę Newtona-Raphsona polegającej na iteracyjnym wyznaczeniu przybliżonej wartości pierwiastka równania,
W(x)=a₀ + a₁ x + a₂ x² + ... + a_n xⁿ
gdzie a₀, a₁, a₂, … a_n należą do zbioru liczb rzeczywistych R i przynajmniej a_n≠0.
Wyznaczenie pierwiastka wymaga określenia punktu startowego x₀ (wartości niewiadomej, od której w sposób iteracyjny będzie poszukiwany pierwiastek). W każdej iteracji dla x_k (aktualnej wartości niewiadomej) obliczana jest wartość wielomianu W(x_k) i w przypadku, gdy W(x_k)≠0 (równanie wielomianu nie jest spełnione), dokonuje się korekty jej wartości. Ideę pokazuje rysunek 1.

proot01.png

Początkowa wartość niewiadomej to x₀, dla której obliczana jest wartość wielomianu. Jeżeli wartość wielomianu W(x_k)≠0, to wystawiana jest styczna do krzywej wielomianu w punkcie x₀. By znaleźć równanie stycznej niezbędna jest wartość pochodnej w tym punkcie W '(x₀). Równanie stycznej to y = W '( x₀ ) x + b, gdzie parametr b jest wyliczany ze współrzędnych punktu styku W( x₀ )=W'( x₀ )x + b → b=W( x₀ )-W '( x₀ )x. Styczna przecina oś X w kolejnym punkcie, którego położenie x₁ można wyznaczyć z równania stycznej: y=0 → 0=W '(x₀)x₁+W(x₀)-W '(x₀) x₀ czyli x₁=x₀ - W(x₀) / W '(x₀. Wyznaczona wartość staje się kolejnym przybliżeniem do rozwiązania.
Dokładnie tą ideę można przenieść do przestrzeni liczb zespolonych, co sprowadza się do rozszerzenia znaczenia argumentu z liczb rzeczywistych do liczb zespolonych. Równocześnie zmienia się wartość funkcji z wartości rzeczywistych na wartości zespolone. Rozwiązaniem jest znalezienie takiego argumentu, gdzie wartość bezwzględna z wartości wielomianu dla tego argumentu jest zerowa.
Jeżeli w trójwymiarowym układzie współrzędnych przedstawić bezwzględną wartości wielomianu w funkcji części rzeczywistej oraz urojonej argumentów zespolonych, to przykładowo wielomian W(x)=x²+1 pokazuje rysunek 2.

proot02.png

Miejsca zerowe wielomianu (jego rozwiązania), to miejsca styku powierzchni wielomianu z płaszczyzną utworzoną przez osie Re z i Im z.
W rzeczywistości w obliczeniach numerycznych za wartość zerową przyjmuje się wartość dostatecznie bliską zera (-ξ < x < ξ ) →x=0). Algorytm Newtona-Raphsona znajduje jeden pierwiastek (jakiś), więc by określić kolejny koniecznością jest „pozbycie się” jego z równania. W tym celu wielomian jest dzielony przez dwumian i obniżany jest jego stopień. W myśl twierdzenia Bézouta wielomian ma obowiązek podzielić się bez reszty.
W przypadku gdy znaleziony pierwiastek x_r jest rzeczywisty, wielomian jest dzielony przez dwumian (x-x_r). W przypadku, gdy rozwiązaniem jest pierwiastek zespolony x_r=x_re + x_im j, to przy istniejących założeniach dotyczących wartości współczynników przy kolejnych potęgach wielomianu, musi istnieć pierwiastek sprzężony x_r=x_re - x_im j. Wielomian jest dzielony przez trójmian x²-2 x_re x + x_re² + x_im² i jego stopień jest obniżony o dwa.
Po obniżeniu stopnia wielomianu ponownie jest poszukiwany kolejny pierwiastek wielomianu aż do osiągnięcia wielomianu drugiego stopnia lub pierwszego stopnia. W tych przypadkach wielomian jest rozwiązywany w oparciu o metody analityczne.
Istotnym elementem w całym zagadnieniu jest optymalne przyjęcie pierwszego przybliżenia rozwiązania (wybranie x₀). Klasyczna metoda Newtona-Raphsona wymaga istnienia rozwiązania rzeczywistego i słabo sobie radzi z rozwiązaniami jedynie zespolonymi (rysunek 3). Uogólniona metoda do przestrzeni zespolonej dobrze zdaje egzamin w takich przypadkach (nie ma wstępnego założenia co do ilości pierwiastków czysto rzeczywistych lub zespolonych). Można przyjąć za pierwsze przybliżenie jednostkę urojoną x=j, jednak czysto urojona wartość początkowa z kolei słabo sobie radzi w sytuacjach, gdzie rozwiązania są czysto rzeczywiste (rysunek 4).

proot03.png

proot04.png

W takich przypadkach program wpada w oscylacje rozwiązania, gdzie kolejne przybliżenia oscylują wzdłuż osi rzeczywistej lub urojonej nie chcąc przejść na drugą oś.
Drogą różnych eksperymentów przyjęcie jako wartości startowej pierwszego przybliżenia liczby x₀=1 + j sprawia, że zauważone problemy znikają.

Rozwiązania w programie

Algorytm poszukiwania miejsc zerowych metodą Newtona-Raphsona został zaimplementowany w języku C w oparciu o kompilator języka C/C++: DEV C++ version 5.11 (narzędzie dostępne na licencji GNU).
Na potrzeby określenia wielomianu oraz jego rozwiązania, w programie został zdefiniowany typ zmiennej.

Typ PolynomialCoeffArrayT jest tablicą współczynników i reprezentuje wielomian (określa współczynniki przy kolejnych potęgach w wielomianie) z tym, że pod indeksem 0 zawiera się wyraz wolny oraz pod indeksem n występuje współczynnik przy n-tej potędze (rysunek 5). Typ PolyRootsArrayT dotyczy rozwiązań wielomianu i jest tablicą o wyrazach zespolonych.

proot09.png

Kompletne rozwiązanie posiłkuje się funkcjami i procedurami pomocniczymi, do których należą przykładowo procedury dzielenia wielomianu.
Do dzielenia wielomianu przez dwumian używana jest procedura DividePoly1.

Dzieli ona wejściowy wielomian określony przez tablicę PolyCoeff stopnia określonego przez parametr PolyDegree przez dwumian (x + PolyRoot). Wynik dzielenia znajduje się w tablicy PolyCoeff.
W przypadku dzielenia wielomianu przez trójmian, używana jest procedura DividePoly2.

Dzieli ona wejściowy wielomian określony przez tablicę PolyCoeff stopnia określonego przez parametr PolyDegree przez trójmian (x² + PolyRootX x + PolyRoot). Wynik dzielenia znajduje się w tablicy PolyCoeff.
Funkcja do poszukiwania pierwiastka jest następująca.

Dostaje on w parametrach dane opisujące wielomian (tablica współczynników PolyCoeff) oraz informacje dotyczące stopnia wielomianu (parametr PolyDegree). Wartość pierwiastka jest wynoszona na zewnątrz poprzez parametr PolyRoot (typu zespolonego). W przypadku pierwiastka rzeczywistego, zawiera się on w części rzeczywistej zmiennej, w przypadku rozwiązania zespolonego, w obu częściach (rzeczywistej oraz urojonej).
Wynikiem funkcji jest status operacji, jako jeden z następujących wariantów.

Ich znaczenie jest następujące:

• LocateRoots_OK – rozwiązanie jest kompletne (wynik funkcji PolyRootsNewtonRaphson),
• LocateRoots_real – funkcja LocateRootsNewtonRaphson doiterowała się pierwiastka rzeczywistego,
• LocateRoots_complex – funkcja LocateRootsNewtonRaphson doiterowała się pierwiastka zespolonego, podany jest jedynie jeden pierwiastek i należy sobie automatycznie dodać pierwiastek sprzężony,
• LocateRoots_fail – wystąpił błąd wewnętrzny (nie powinno się zdarzyć),
• LocateRoots_overflow – w poszukiwaniu pierwiastka doszło do dzielenia przez zero (normalnie nie powinno wystąpić),
• LocateRoots_LoopLimit - funkcja LocateRootsNewtonRaphson wpadła w oscylacje i nie doiterowała się żadnego rozwiązania (limit iteracji jest określony przez stałą LoopLimit.

Rozwiązanie funkcji poszukiwania pierwiastka jest następujące.

W powyższej funkcji używana jest pomocnicza funkcja do obliczenia wartości wielomianu oraz wartości pochodnej dla określonego argumentu. Jej postać pokazuje:

Idea obliczeń jest oparta o następujące operacje realizowane w pętli:

• FValue jest podstawiona na wartość współczynnika przy najwyższej potędze a_n (FV=a_n),
  wykonane są pętli następujące operacje (arg jest argumentem funkcji → parametr wywołania),
  FV= a_n arg + a_n-1
  FV=( a_n arg + a_n-1) arg + a_n-2
  …
  FV=((( a_n arg + a_n-1) arg + a_n-2) arg + a_n-3 ... ) arg + a₀

Do obliczania wartości pochodnej używana jest również powyższa funkcja, jej realizację pokazuje:

Finalnie do znalezienia wszystkich pierwiastków wielomianu używana jest funkcja PolyRootsNewtonRaphson .

Jej zadaniem jest obniżanie stopnia wielomianu aż do uzyskania końcowego rozwiązania.
Przykład użycia pokazuje:


Dla tropicieli:

Testy

Wielomian trzeciego stopnia W(x)=3 x³ + 3 x² + 3 x. Z samej matematyki wiadomo, że wielomian ma pierwiastek rzeczywisty o wartości x_r=0, gdzie po podzieleniu wielomianu przez dwumian powstaje równanie kwadratowe, które rozwiązane jest już analitycznie. Znalezienie pierwszego pierwiastka odbywa się iteracyjnie. Program wykonał następujące korki (rysunek 6).

proot05.png

Wielomian W(x)=x⁴ – 1000256 x² + 256000000=(x² - 1000000)(x²-256), który ma cztery pierwiastki rzeczywiste. Poszukując pierwszego pierwiastka program wykonał następujące iteracje (rysunek 7).

proot06.png

Poszukiwanie drugiego pierwiastka równania generuje następujące działania (rysunek 8). Po znalezieniu tego pierwiastka, równanie staje się trójmianem i dalej jest rozwiązywane analitycznie.

proot07.png

Bardziej złożony test dotyczy następującego wielomianu:
W ( x ) = x¹⁹ – x¹⁸ – 2 x¹⁷ -3 x¹⁶ – 4 x¹⁵ -1 x¹⁴ + x¹³ + x¹² – x¹⁰ – 2 x⁹ -3 x⁸ – x⁷ -5 x⁶ +5 x⁵ + 4 x⁴ + 3x³
Zbieżność poszczególnych rozwiązań pokazuje rysunek 9 (w przypadku rozwiązań zespolonych, rysunek nie zawiera drogi dojścia dla obu pierwiastków sprzężonych: tylko jednego z nich). Wykonane cykle iteracyjne dla poszczególnych pierwiastków prezentuje listing.

proot08.png

Ostatni pierwiastek wyliczony jest metodą analityczną. Finalnie równanie ma następujący zestaw pierwiastków:

Generator liczb losowych

W wielu zagadnieniach technicznych przydatnym jest generator liczb losowych. Za jego pomocą można przykładowo symulować sygnał zakłócający i badać metody jego oczyszczania oraz opracowywać metody eliminacji. Jedną z istotnych cech generatora liczb losowych jest jego kolor. To dosyć niepoważnie brzmiące określenie ma głęboki sens fizyczny. Wyprodukowany sygnał można poddać analizie częstotliwościowej. Jeżeli zawartość harmonicznych będzie dosyć równomiernie rozłożona, to można mówić o białym generatorze szumu (zawiera każdą możliwą częstotliwość). Jeżeli tak nie jest, to generator produkuje przykładowo kolor różowy. Uzyskanie generatora o określonym kolorze nie należy do czynności trywialnych. Jest sporo wręcz sprzecznych ze sobą wymagań, jak choćby nieskończenie wielki okres (z tą nieskończonością to chyba lekkie przegięcie), czyli po ilu losowaniach powtórzą się wartości. W literaturze można znaleźć wiele przykładowych generatorów i uważam, że nie warto traciś sił i środków na wymyślanie własnego. Poniższy jest zaczerpnięty z "Numerical Recipes in C". Typowe generatory dają wartości jako liczby zmiennoprzecinkowe z zakresu 0 .. 1.

Uruchomienie tego programu generuje plik tekstowy z generowanymi wartościami, jego fragment (lektura całości jest raczej monotonna )

Numeryczne rozwiązywanie liniowych
układów równań

W wielu zagadnieniach technicznych czynnością wykonywaną niejako w tle, jest rozwiązanie układu równań liniowych. To element składowy wielu działań, taki mały trybik w większej maszynerii. Równania można zapisać w klasycznej formie:

mat01_form00.png

lub w formie macierzowej:

mat01_form01.png

gdzie A jest macierzą kwadratową n*n układu składającą się ze wszystkich współczynników układu równań,
X jest poszukiwaną macierzą kolumnową
B jest macierzą kolumnową współczynników (tych po drugiej stronie równań).

mat01_form02.png

mat01_form03.png

mat01_form04.png

By rozwiązać ten układ równań, należy go pomnożyć obustronnie przez macierz odwrotną do macierzy A (pamiętając, że działania należy wykonać we właściwej kolejności, gdyż nie są przemienne).

mat01_form05.png

Ta prosta koncepcja zawiera w sobie kilka mniejszych, że wymienię przykładowo mnożenie macierzy kwadratowej przez macierz kolumnową czy odwracanie macierzy. Tu z kolei elementem podstawowym jest obliczenie wyznacznika macierzy. Zagłębiając się w cały ten ciąg operacji podstawowym elementem jest obliczenie wyznacznika. W przypadku macierzy 2*2, wynik jest nawet prosty do zapamiętania:

mat01_form06.png

w przypadku 3*3, to już staje się bardziej złożone (wiadomo, im bardziej w las, tym więcej drzew), ale jeszcze do opanowania:

mat01_form07.png

W bardziej złożonych sytuacjach to zaczynają się już „schody”. Jednak tu znanych jest kilka możliwości i naszym zadaniem jest odpowiedni wybór: spośród iluś możliwości wybieramy jedną. Zdecydowałem się na metodę Laplace'a, gdyż jest samopodobną (rozkłada się na ciąg działań o identycznej budowie jednak obejmujących mniejsze struktury).
Przykładowo dla macierzy 4*4 rozkłada się na obliczenie czterech wyznaczników macierzy 3*3 powstałych poprzez wykreślenie jednej kolumny i odpowiednich wierszy. Ja wybrałem w implementacji programowej pierwszą kolumnę i odpowiednie wiersze (rys 1).

matrix00.png

W przypadku ręcznego rachowania to można dokonywać pewnych optymalizacji (przykładowo wybrać takie wiersze/kolumny, gdzie występuje dużo zer). Z punktu widzenia programu to ma mniejsze znaczenie, gdyż wtedy koniecznością staje się stworzenie bardziej złożonych algorytmów, które pozwolą na jakąś tam optymalizację (coś jak sztuczna inteligencja?). Narobić się w programie, by praktycznie nic nie uzyskać, to może nie jest najlepsza koncepcja. Dla procka to jedna rybka, co robi, umie mnożyć, to niech mnoży (robi to wystarczająco dobrze). No więc implementacja jest bez wodotrysków.
Oczywiście obliczenie wyznaczników o większych rozmiarach sprowadza się do obliczenie iluś tam wyznaczników macierzy mniejszych, które z kolei rekurencyjnie składają się z obliczenia kolejnych wyznaczników macierzy mniejszych itd.

Program

Utworzona jest następująca struktura do „trzymania” macierzy:

i funkcja obliczająca wartość wyznacznika:

Przykład użycia i testy:


Puszczenie programu dla takich danych generuje następujący rezultat:


Dla tropicieli:

Numeryczne rozwiązywanie liniowych
układów równań - cd

Rozwiązanie układu równań wymaga znalezienia macierzy odwrotnej. Jest to jakaś tam złożona operacja, gdzie również istnieją dwie drogi postępowania: metoda dopełnień algebraicznych (w ogromnej mierze opartej na obliczaniu wyznaczników) oraz metoda eliminacji Gaussa. Ponieważ mamy już stworzone narzędzia (liczenie wyznaczników), implementacja będzie opierała się o dopełnienia algebraiczne.
Dopełnienie określonego elementu macierzy to wyznacznik macierzy powstałej z wykreślenia z macierzy wiersza i kolumny wziętej z odpowiednim znakiem (wynikającym z parzystości/nieparzystości położenia tego elementu).
Nie będę rozwodził się na temat szczegółów metody, gdyż ta wiedza jest powszechnie dostępna i można samodzielnie z nią się zapoznać. Pierwszym krokiem, jaki należy wykonać, to obliczyć wyznacznik macierzy. Może się tak zdarzyć, że wyznacznik macierzy wyjdzie zero, to wtedy jest już pozamiatane: nie jest możliwe obliczenie macierzy odwrotnej a sam układ równań jest nierozwiązywalny.
Opierając się o wyżej opisaną funkcję obliczania wyznaczników, implementacja algorytmu znajdowania macierzy odwrotnej jest następująca:

gdzie pomocnicze elementy są następujące (ustalenie znaku dopełnienia):

dopełnienie algebraiczne macierzy (usunięcie z macierzy niepożądanych elementów):

oraz transponowanie macierzy (zamiana elementów miejscami):

Do pełni szczęścia konieczna jest jeszcze operacja mnożenia macierzy kwadratowej przez macierz kolumnową. Dostarczam gotowego rozwiązania:

Pozostało sprawdzić co z tego wyszło. Kawałek programu testowego:

W konkretnym przypadku rozwiązywany jest układ równań:

• -1 x₁ + 2 x₂ + 3 x₃ + 4 x₄ + 5 x₅ + 6 x₆ = 19
• 11 x₁ - 12 x₂ + 13 x₃ + 14 x₄ + 15 x₅ + 16 x₆ = 57
• 21 x₁ + 22 x₂ - 23 x₃ + 24 x₄ + 25 x₅ + 26 x₆ = 95
• 31 x₁ + 32 x₂ + 33 x₃ - 34 x₄ + 35 x₅ + 36 x₆ = 133
• 41 x₁ + 42 x₂ + 43 x₃ + 44 x₄ - 45 x₅ + 46 x₆ = 171
• 51 x₁ + 52 x₂ + 53 x₃ + 54 x₄ + 55 x₅ - 56 x₆ = 209

Zostaje uruchomiony program, który wyprodukował następujący wynik:


Dla tropicieli: