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)=a0 + a1 x + a2 x2 + ... + an xn
gdzie a0, a1, a2, … an należą do zbioru liczb rzeczywistych R i przynajmniej an≠0.
Wyznaczenie pierwiastka wymaga określenia punktu startowego x0 (wartości niewiadomej, od której w sposób iteracyjny będzie poszukiwany pierwiastek). W każdej iteracji dla xk (aktualnej wartości niewiadomej) obliczana jest wartość wielomianu W(xk) i w przypadku, gdy W(xk)≠0 (równanie wielomianu nie jest spełnione), dokonuje się korekty jej wartości. Ideę pokazuje rysunek 1.
Początkowa wartość niewiadomej to x0, dla której obliczana jest wartość wielomianu. Jeżeli wartość wielomianu W(xk)≠0, to wystawiana jest styczna do krzywej wielomianu w punkcie x0. By znaleźć równanie stycznej niezbędna jest wartość pochodnej w tym punkcie W '(x0). Równanie stycznej to y = W '( x0 ) x + b, gdzie parametr b jest wyliczany ze współrzędnych punktu styku W( x0 )=W'( x0 )x + b → b=W( x0 )-W '( x0 )x. Styczna przecina oś X w kolejnym punkcie, którego położenie x1 można wyznaczyć z równania stycznej: y=0 → 0=W '(x0)x1+W(x0)-W '(x0) x0 czyli x1=x0 - W(x0) / W '(x0. 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)=x2+1 pokazuje rysunek 2.
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 xr jest rzeczywisty, wielomian jest dzielony przez dwumian (x-xr). W przypadku, gdy rozwiązaniem jest pierwiastek zespolony xr=xre + xim 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 xr=xre - xim j. Wielomian jest dzielony przez trójmian x2-2 xre x + xre2 + xim2 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 x0). 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).
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 x0=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.
Kod: Zaznacz cały
#define MaxPolynomialDegree 20
typedef double PolynomialCoeffArrayT [ MaxPolynomialDegree ] ;
// PolynomialCoeffArrayT [ 0 ] - x^0
// PolynomialCoeffArrayT [ 1 ] - x^1
// PolynomialCoeffArrayT [ 2 ] - x^2
// .....
// PolynomialCoeffArrayT [ n ] - x^n
typedef double complex PolyRootsArrayT [ MaxPolynomialDegree - 1 ] ;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.
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.
Kod: Zaznacz cały
void DividePoly1 ( double PolyCoeff [ ] ,
double PolyRoot ,
unsigned short PolyDegree )
{
signed short Loop ;
PolynomialCoeffArrayT QuotientPoly ;
/*------------------------------------------------------------------*/
for ( Loop = PolyDegree - 1 ; Loop >= 0 ; Loop -- )
{
QuotientPoly [ Loop ] = PolyCoeff [ Loop + 1 ] ;
PolyCoeff [ Loop ] = PolyCoeff [ Loop ] - PolyRoot * QuotientPoly [ Loop ];
} /* for */ ;
for ( Loop = 0 ; Loop < PolyDegree ; Loop ++ )
PolyCoeff [ Loop ] = QuotientPoly [ Loop ] ;
} /* 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.
Kod: Zaznacz cały
void DividePoly2 ( double PolyCoeff [ ] ,
double PolyRootX ,
double PolyRoot ,
unsigned short PolyDegree )
{
signed short Loop ;
PolynomialCoeffArrayT QuotientPoly ;
/*------------------------------------------------------------------*/
for ( Loop = PolyDegree - 2 ; Loop >= 0 ; Loop -- )
{
QuotientPoly [ Loop ] = PolyCoeff [ Loop + 2 ] ;
PolyCoeff [ Loop + 1 ] = PolyCoeff [ Loop + 1 ] - PolyRootX * QuotientPoly [ Loop ] ;
PolyCoeff [ Loop ] = PolyCoeff [ Loop ] - PolyRoot * QuotientPoly [ Loop ] ;
} /* for */ ;
for ( Loop = 0 ; Loop < PolyDegree ; Loop ++ )
PolyCoeff [ Loop ] = QuotientPoly [ Loop ] ;
} /* DividePoly2 */Dzieli ona wejściowy wielomian określony przez tablicę PolyCoeff stopnia określonego przez parametr PolyDegree przez trójmian (x2 + PolyRootX x + PolyRoot). Wynik dzielenia znajduje się w tablicy PolyCoeff.
Funkcja do poszukiwania pierwiastka jest następująca.
Kod: Zaznacz cały
unsigned short LocateRootsNewtonRaphson ( double complex * PolyRoot ,
double PolyCoeff [ ] ,
unsigned short PolyDegree )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.
Kod: Zaznacz cały
#define LocateRoots_OK 0
#define LocateRoots_real 1
#define LocateRoots_complex 2
#define LocateRoots_fail 3
#define LocateRoots_overflow 4
#define LocateRoots_LoopLimit 5Ich 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.
Kod: Zaznacz cały
unsigned short LocateRootsNewtonRaphson ( double complex * PolyRoot ,
double PolyCoeff [ ] ,
unsigned short PolyDegree )
{
double Delta ;
double complex Root ;
double complex CurrArg ;
double complex PolyValue ;
double complex DerivValue ;
unsigned short LoopCounter ;
/*------------------------------------------------------------------*/
if ( PolyDegree < 1 )
return ( LocateRoots_fail ) ;
if ( PolyDegree == 1 )
{
if ( IsZero ( PolyCoeff [ 1 ] ) )
return ( LocateRoots_overflow ) ;
__real__ Root = - PolyCoeff [ 0 ] / PolyCoeff [ 1 ] ;
__imag__ Root = 0.0 ;
* PolyRoot = Root ;
return ( LocateRoots_real ) ;
} /* if */ ;
if ( PolyDegree == 2 )
{
if ( IsZero ( PolyCoeff [ 2 ] ) )
return ( LocateRoots_overflow ) ;
Delta = PolyCoeff [ 1 ] * PolyCoeff[1] - 4.0 * PolyCoeff [ 2 ] * PolyCoeff [0] ;
if ( Delta >= 0 )
{
__real__ Root = ( -PolyCoeff[1] + sqrt ( Delta )) / ( 2.0 * PolyCoeff [2] ) ;
__imag__ Root = 0.0 ;
* PolyRoot = Root ;
return ( LocateRoots_real ) ;
} /* if ... */
else
{
__real__ Root = - PolyCoeff [ 1 ] / ( 2.0 * PolyCoeff [ 2 ] ) ;
__imag__ Root = - sqrt ( - Delta ) / ( 2.0 * PolyCoeff [ 2 ] ) ;
* PolyRoot = Root ;
return ( LocateRoots_complex ) ;
} /* if ... else */ ;
return ( LocateRoots_fail ) ;
} /* if */ ;
CurrArg = 1.0 + I ;
LoopCounter = 0 ;
for ( ; ; )
{
PolyValue = FunctionValue ( PolyCoeff , PolyDegree , CurrArg ) ;
if ( IsZero ( __real__ PolyValue ) && IsZero ( __imag__ PolyValue ) )
{
if ( IsZero ( __imag__ CurrArg ) )
{
__imag__ CurrArg = 0.0 ;
* PolyRoot = CurrArg ;
return ( LocateRoots_real ) ;
} /* if ... */
else
{
* PolyRoot = CurrArg ;
return ( LocateRoots_complex ) ;
} /* if ... else */ ;
} /* if */ ;
DerivValue = DerivativeValue ( PolyCoeff , PolyDegree , CurrArg ) ;
Root = CurrArg ;
if ( IsZero ( __real__ DerivValue ) && IsZero ( __imag__ DerivValue ) )
CurrArg = 0.9 * Root ;
else
CurrArg = Root - PolyValue / DerivValue ;
LoopCounter ++ ;
if ( LoopCounter > LoopLimit )
return LocateRoots_LoopLimit ;
} /* for */ ;
} /* LocateRootsNewtonRaphson */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:
Kod: Zaznacz cały
double complex FunctionValue ( double PolyCoeff [ ] ,
unsigned short PolyDegree ,
double complex Argument )
{
double complex FValue ;
signed short Loop ;
/*------------------------------------------------------------------*/
FValue = PolyCoeff [ PolyDegree ] + 0.0 * I ;
for ( Loop = PolyDegree ; Loop > 0 ; Loop -- )
{
FValue = FValue * Argument ;
__real__ FValue = __real__ FValue + PolyCoeff [ Loop - 1 ] ;
} /* for */ ;
return ( FValue ) ;
} /* FunctionValue */ 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 an (FV=an),
wykonane są pętli następujące operacje (arg jest argumentem funkcji → parametr wywołania),
FV= an arg + an-1
FV=( an arg + an-1) arg + an-2
…
FV=((( an arg + an-1) arg + an-2) arg + an-3 ... ) arg + a0
Kod: Zaznacz cały
double complex DerivativeValue ( double PolyCoeff [ ] ,
unsigned short PolyDegree ,
double complex Argument )
{
signed short Loop ;
PolynomialCoeffArrayT Derivative ;
/*------------------------------------------------------------------*/
for ( Loop = PolyDegree ; Loop > 0 ; Loop -- )
Derivative [ Loop - 1 ] = ( double ) ( Loop ) * PolyCoeff [ Loop ] ;
return ( FunctionValue ( Derivative , PolyDegree - 1 , Argument ) ) ;
} /* DerivativeValue */Finalnie do znalezienia wszystkich pierwiastków wielomianu używana jest funkcja PolyRootsNewtonRaphson .
Kod: Zaznacz cały
unsigned short PolyRootsNewtonRaphson ( double complex PolyRootArr [ ] ,
double PolyCoeff [ ] ,
unsigned short PolyDegree )
{
double complex PolyRoot ;
unsigned short RootsIndex ;
/*------------------------------------------------------------------*/
RootsIndex = 0 ;
for ( ; ; )
{
switch ( LocateRootsNewtonRaphson ( & PolyRoot , PolyCoeff , PolyDegree ) )
{
case LocateRoots_real :
DividePoly1 ( PolyCoeff , - __real__ PolyRoot , PolyDegree ) ;
PolyDegree -- ;
PolyRootArr [ RootsIndex ++ ] = PolyRoot ;
break ;
case LocateRoots_complex :
DividePoly2 ( PolyCoeff , - 2.0 * __real__ PolyRoot ,
__real__ PolyRoot * __real__ PolyRoot +
__imag__ PolyRoot * __imag__ PolyRoot , PolyDegree ) ;
PolyDegree -- ;
PolyDegree -- ;
PolyRootArr [ RootsIndex ++ ] = PolyRoot ;
__imag__ PolyRoot = - __imag__ PolyRoot ;
PolyRootArr [ RootsIndex ++ ] = PolyRoot ;
break ;
case LocateRoots_fail :
return ( LocateRoots_fail ) ;
case LocateRoots_overflow :
return ( LocateRoots_overflow ) ;
case LocateRoots_LoopLimit :
return ( LocateRoots_LoopLimit ) ;
} /* switch */ ;
if ( PolyDegree == 0 )
return ( LocateRoots_OK ) ;
} /* for */ ;
} /* PolyRootsNewtonRaphson */Jej zadaniem jest obniżanie stopnia wielomianu aż do uzyskania końcowego rozwiązania.
Przykład użycia pokazuje:
Kod: Zaznacz cały
int main ( int argc ,
char * argv [ ] )
{
unsigned short Loop ;
unsigned short PDegree ;
double complex FValue ;
PolynomialCoeffArrayT PolyCopy ;
/*------------------------------------------------------------------*/
OutFile = fopen ( "proots.txt" , "w" ) ;
fprintf ( OutFile , "*** Program poszukujacy miejsc zerowych wielomianu ***\n\n" ) ;
Polynomial [ 0 ] = 0.0 ;
Polynomial [ 1 ] = 3.0 ;
Polynomial [ 2 ] = 3.0 ;
Polynomial [ 3 ] = 3.0 ;
PDegree = 3 ;
for ( Loop = 0 ; Loop <= PDegree ; Loop ++ )
PolyCopy [ Loop ] = Polynomial [ Loop ] ;
PrintPoly ( Polynomial , PDegree ) ;
PolyRootsNewtonRaphson ( PolyRoots , Polynomial , PDegree ) ;
fprintf ( OutFile , "\n\n\nWyniki:\n" );
for ( Loop = 0 ; Loop < PDegree ; Loop ++ )
{
if ( IsZero ( __imag__ PolyRoots [ Loop ] ) )
fprintf ( OutFile , "Pierw.rzecz.:(%.10f)" , __real__ PolyRoots [ Loop ] ) ;
else
fprintf ( OutFile , "Pierw.zespo.: (%.10f)+(%.10f j)" ,
__real__ PolyRoots [ Loop ] , __imag__ PolyRoots [ Loop ] ) ;
FValue = FunctionValue ( PolyCopy , PDegree , PolyRoots [ Loop ] ) ;
fprintf ( OutFile , " wart.wielom.: (%.10f)+(%.10f j)\n" ,
__real__ FValue , __imag__ FValue ) ;
} /* for */ ;
fclose ( OutFile ) ;
return 0 ;
} /* main */
Dla tropicieli:
)