Wielomiany w octavie
W octavie istnieje cała gamma funkcji związanych z wielomianami:
-
polyval( )
- funkcja pozwalająca obliczać wartość wielomianu
zadanego w bazie potęgowej dla danej wartości, czy całej macierzy wartości
-
polyfit( )
- funkcja znajdująca współczynniki wielomianu zadanego stopnia najlepiej dopasowanego do zadanej tabelki punktów. -
polyinteg( )
- zwraca współczynniki wielomianu będącego całką nieoznaczoną z danego wielomianu -
polyderiv( )
- zwraca współczynniki wielomianu będącego pochodną z danego wielomianu -
roots( )
- zwraca zera wielomianu -
conv(a,b)
- zwraca współczynniki wielomianu będącego iloczynem wielomianów o współczynnikach odpowiednio z wektorów
Będziemy korzystali przede wszystkim z dwóch pierwszych funkcji.
We wszystkich tych funkcjach przyjmowane jest, że rozpatrujemy współczynniki wielomianu w bazie potęgowej, ale w następującej kolejności:
Współczynniki indeksowane są od jedynki, tzn. wielomian stopnia nie większego od :
jest reprezentowany przez wektor współczynników: .
Jeśli chcemy policzyć wartość w danym punkcie
, czy ewentualnie dla tablicy punktów umieszczonych w macierzy
,
wywołujemy polyval(a,x)
lub polyval(a,X)
.
Tak więc, aby narysować wykres funkcji na
można wykonać następującą sekwencję komend :
Interpolacja wielomianowa
Interpolacja wielomianowa polega na tym, że szukamy funkcji z pewnej przestrzeni wielomianów (zazwyczaj wielomianów stopnia nie większego od
), które w tych punktach spełniają odpowiednie warunki interpolacyjne, tzn.
oraz jej pochodne mają w tych punktach zadane wartości.
Interpolacja Lagrange'a
W interpolacji wielomianowej Lagrange'a dla zadanych różnych punktów
i wartości
szukamy wielomianu
stopnia co najwyżej
takiego, że spełnione są następujące warunki interpolacyjne:
Funkcja polyfit(x,y,n)
znajduje
współczynniki wielomianu w bazie potęgowej dla wektorów
długości
. Ważne aby wartości w
były różne.
Przetestujmy tę funkcję dla następujących danych:
,
:
Widać na rysunku [link], że wielomian spełnia warunki interpolacyjne.
Możemy to sprawdzić też obliczeniowo:
Otrzymaliśmy wynik: błąd w węzłach jest równy prawie zero.
Przetestujmy polyfit(x,y,n)
dla dużych dla węzłów równoodległych na
i losowe wartości
z
, następnie policzmy błąd w węzłach w zależności od
:
Widzimy, że możemy otrzymać niepoprawny wynik dla dużych .
Wynika to ze złych własności algorytmu stosowanego przez octave'a ze względu na zaburzenia powodowane przez niedokładne obliczenia w arytmetyce zmiennopozycyjnej.
Zbadajmy błąd pomiędzy funkcją, a jej wielomianem interpolacyjnym.
Na początek rozpatrzmy analityczną funkcję i węzły równoodległe, oraz normę supremum na odcinku
. Normę supremum:
przybliżymy poprzez dyskretną normę na siatce równomiernie rozłożonych punktów na tym odcinku
, czyli
z :
Interpolacja Hermite'a
Rozpatrzmy teraz zadanie interpolacji Hermite'a.
W interpolacji wielomianowej Hermite'a dla zadanych różnych punktów
i naturalnych krotności węzłów
szukamy wielomianu
stopnia co najwyżej
dla
takiego, że spełnione są następujące warunki interpolacyjne:
Tutaj to
zadanych wartości.
W octavie nie ma funkcji realizującej interpolacje Hermite'a.
Ale czy w szczególnych przypadkach nie możemy łatwo rozwiązać tego problemu korzystając z odpowiednich funkcji octave'a?
Rozpatrzmy przypadek węzłów tej samej krotności; np. różnych dwukrotnych węzłów.
Chcemy znaleźć współczynniki wielomianu dla
takie, że
Czyli rozwiązanie spełnia następujący układ równań liniowych:
\begin{eqnarray*}
\sum_{k=0}^N a_kx_j^k=f(x_k) && j=0,\ldots,n \\
\sum_{k=1}^N ka_kx_j^{k-1}=f'(x_k) && j=0,\ldots,n.
\end{eqnarray*}
Utwórzmy macierz tego układu z pomocą funkcji vander(x)
, która tworzy macierz Vandermonde'a dla węzłów podanych w wektorze :
#tworzy wspolczynniki wielomianu Hermite'a dla
#wezlow dwukrotnych z x
#y= wartosci w wezlach x (wektor pionowy)
#wartosci pochodnej w wezlach x (wektor pionowy)
#output: wspolczynniki wielomainu w bazie potegowej
# (Zgodne z polyval())
#C- opcja - macierz ukladu
n=length(x);
N=2*n-1;
A=vander(x,N+1);
D=diag(N:-1:0);
B=zeros(size(A));
B(:,1:N)=A(:,2:N+1);
B=B*D;
C=[A;B];
f=[y;dy];
a=C\[y;dy];
endfunction
Na początku przetestujmy tę funkcję
dla węzłów dwukrotnych i funkcji
.
Z wykresu funkcji i wielomianu widzimy, że
wielomian przecina wykres i jest styczny w punktach , por. rysunek [link].
Zobaczmy co się stanie w przypadku trzech węzłów.
Testujemy funkcję dla węzłów dwukrotnych i
.
Z rysunku [link] widać, że funkcja działa poprawnie również w tym przypadku.
%%
% SECTION SPLAJNY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
Interpolacja zespolona. Algorytm FFT
(#)
%%
% SECTION FFT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W tym rozdziale omówimy krótko interpolację zespoloną, ale tylko w przypadku
określonych węzłów zespolonych równomiernie rozmieszczonych na sferze jednostkowej.
Chcemy znaleźć współczynniki zespolone ,
takie, że
dla danych zespolonych i
dla
, czyli pierwiastków z jedynki stopnia
.
Równoważnie możemy to zadanie sformułować jako zadanie znalezienia takich, że
Okazuje się, że rozwiązanie możemy wyrazić poprzez operator dyskretnej transformaty Fouriera, czy równoważnie jako mnożenie przez macierz :
gdzie , a
.
Mnożenie przez macierz można szybko wykonać z wykorzystaniem algorytmu szybkiej transformacji Fouriera, czyli FFT (Fast Fourier Transform). Odpowiednia wersja algorytmu FFT służy też szybkiemu mnożeniu przez macierz odwrotną do
:
.
Obliczając wartość wielomianu w punktach
, czy równoważnie wartości wielomianu trygonometrycznego
w punktach
, to przyjmując
otrzymujemy:
Warto dodać, że często macierz DFT definiuje się bez czynnika , oczywiście wtedy macierz odwrotna też musi być odpowiednio przeskalowana.
Funkcja
służy w octave'ie mnożeniu przez macierz
Jej najprostsze wywołanie to
Sprawdźmy, czy funkcja fft()
oblicza wartość DFT zgodnie z naszą definicją. Policzmy , powinniśmy otrzymać
wektor samych jedynek:
a otrzymaliśmy:
4 + 0i
4 + 0i
4 + 0i
4 - 0i
To oznacza, że funkcja fft()
oblicza wartość mnożenia przez
.
Z kolei funkcja ifft()
powinna obliczać mnożenie przez macierz odwrotną do .
Sprawdźmy, jak to działa:
Powtórzmy to dla większego :
Stwórzmy macierz i porównajmy szybkość mnożenia przez tę macierz wykonaną za pomocą standardowego operatora octave'a, czyli operatora
*
, z szybkością działania funkcji fft()
.
Stwórzmy najpierw funkcję tworzącą macierz :
Przetestujmy ją na początek dla :
W tym przypadku wyniki się pokrywają.
Zgodnie z teorią, por. np. [Jankowska:1981:MN], koszt obliczenia DFT przy zastosowaniu algorytmu FFT to
, a standardowe mnożenia przez macierz
kosztuje
.
A teraz testujemy szybkość:
Na moim komputerze FFT było ponad razy szybsze dla
.
Policzmy jeszcze wykres realnego kosztu względem :
Na wynik chwilę musimy poczekać, ale potwierdza on teorię, że FFT jest wyraźnie szybszym algorytmem.
Interpolacja splajnowa
(#)
W przypadku interpolacji splajnowej rozpatrujemy dane węzły:
na odcinku , oraz
funkcje splajnowe (splajny), czyli funkcje, które są odpowiedniej klasy gładkości (czyli są w , oraz obcięte do dowolnego pododcinka
dla
, są wielomianami co najwyżej ustalonego stopnia.
W przypadku splajnów kubicznych rozważamy funkcje, które są klasy na
, oraz są wielomianami kubicznymi na każdym pododcinku.
Zadanie interpolacji splajnowej kubicznej polega na znalezieniu splajnu
kubicznego takiego, że
dla zadanych ,
oraz spełniającego odpowiednie warunki brzegowe. Np. w przypadku splajnu hermitowskiego
musi spełnić warunki brzegowe hermitowskie:
dla zadanych dodatkowych dwóch wartości .
W przypadku splajnu naturalnego warunki brzegowe to
Rozpatruje się również splajny typu not-a-knot. W tym przypadku,
zamiast warunków brzegowych, dodaje się sztucznie dwa warunki ciągłości trzeciej pochodnej w węzłach i
.
Z kolei splajny okresowe spełniają następujące warunki brzegowe okresowe:
W każdym z tych czterech przypadków splajn interpolacyjny jest wyznaczony jednoznacznie, por. np.
[Kincaid:2006:AN].
W octave'ie istnieje kilka funkcji związanych z interpolacją splajnami kubicznymi, czy ogólnie - z funkcjami wielomianowymi na pododcinkach:
-
spline()
- funkcja służąca
znalezieniu splajnu interpolacyjnego hermitowskiego, czy
typu not-a-knot -
ppval()
- funkcja służąca obliczeniu wartości splajnu zadanego w formacie octave'a -
mkpp()
- tworzy funkcję wielomianową na pod-odcinkach (szczegółyhelp mkpp()
) -
unmkpp()
- ze struktury splajnu w octave'ie zwraca współczynniki wielomianów na pod-odcinkach (szczegółyhelp unmkpp()
)
Podstawową funkcją służącą
znalezieniu splajnu interpolacyjnego hermitowskiego, czy
typu not-a-knot jest funkcja spline()
.
Jej najprostsze wywołanie to
gdzie to wektor wymiaru
z węzłami interpolacji splajnowej, a wektor
ma tę samą długość co
i zawiera wartości jakie ma przyjąć splajn w tych węzłach.
Druga możliwość to - przy takim samym wektorze węzłów , podanie wektora
o długości
. Wtedy pierwsza i ostatnia wartość wektora
to wartości pochodnych splajnu w końcowych węzłach
. Pozostałe wartości
tzn.
dla
zawierają wartości splajnu w węzłach z
, czyli są takie same jak w pierwszym przypadku.
Funkcja zwraca strukturę typu , czyli w odpowiednim formacie octave'a splajnu kubicznego typu not-a-knot w pierwszym przypadku, a w drugim - splajnu hermitowskiego. Następnie, korzystając z funkcji
ppval()
można obliczyć wartość tego splajnu w punkcie, czy tablicy punktów.
Policzmy splajn który w węzłach przyjmie wartości
i narysujmy jego wykres, por. rysunek [link]:
Korzystając z tej samej funkcji, stwórzmy splajn hermitowski przyjmując, że jego pochodne w końcach wynoszą zero. Następnie narysujmy wykresy obu splajnów na odcinku , por. rysunek [link], oraz blisko lewego końca, por. rysunek [link]:
Widać, że splajny są różne, oraz że splajn hermitowski ma pochodną równą zero w lewym i prawym końcu.
Czy w octave'ie można wyznaczyć w prosty sposób inne splajny, np. naturalny lub okresowy?
Okazuje się, że tak, ale trzeba wykorzystać funkcję z rozszerzenia octave'a, czyli octave-forge'a: csape()
.
Jej wywołanie to
gdzie to wektory długości
z węzłami i wartościami splajnu w węzłach, a cond przyjmuje wartości:
- 'variational' w przypadku splajnu interpolacyjnego kubicznego naturalnego
- 'complete' w przypadku splajnu interpolacyjnego kubicznego hermitowskiego,
wartości pochodnych w końcach są podane w parametrzevalc
- 'not-a-knot' w przypadku splajnu interpolacyjnego kubicznego typu \textsl{not-a-knot}
- 'periodic'
w przypadku splajnu interpolacyjnego kubicznego okresowego - 'second'
w przypadku splajnu interpolacyjnego kubicznego z ustalonymi wartościami drugiej pochodnej w końcach, zgodnymi z tym, co podano w parametrzevalc
Funkcja zwraca strukturę typu , czyli dane splajnu kubicznego w formacie octave'a.
Porównajmy wyniki tej funkcji z wynikiem funkcji spline()
dla danych z naszego prostego przykładu i splajnu typu not-a-knot. Policzymy dyskretną normę maksimum (na siatce równomiernej 300 punktów na z wyników obu funkcji:
Otrzymaliśmy zero.
Na jednym wykresie narysujmy wykresy trzech splajnów interpolacyjnych dla z tymi samymi warunkami interpolacyjnymi, ale z różnymi warunkami brzegowymi: hermitowskim, not-a-knot i naturalnym.
Wszystkie trzy splajny są różne w tym przypadku, por.
rysunek [link].
Oczywiście obie funkcje umożliwiają dowolny wybór węzłów, np. weźmy węzły
z tymi samymi wartościami i stwórzmy splajny obu typów, por. rysunek [link]:
Popatrzmy na błędy aproksymacji w normie supremum.
Ustalmy, że interpolujemy splajnami kubicznymi znaną funkcję gładką na czterech węzłach równoodległych na odcinku
.
Jeśli to struktura typu
opisująca splajn interpolujący
, to dyskretną normę maksimum różnicy między
, a splajnem
, np. na siatce o tysiącu punktach, możemy obliczyć komendą:
Założyliśmy, że funkcja jest zaimplementowana wektorowo, tzn. że wywołanie w octave'ie
f(z)
dla wektora zwróci wektor
. W przeciwnym razie należałoby użyć pętli do wyznaczenia wektora
:
W poniższym kodzie obliczymy przybliżoną normę maksimum na odcinku pomiędzy
, a jej dwoma splajnami interpolacyjnymi kubicznymi: naturalnym i typu not-a-knot na czterech węzłach równoodległych:
Błąd dla splajnu naturalnego był większy:
errn = 0.048190
, niż dla splajnu typu not-a -knot:
err = 0.025120
, co widać też na rysunku [link].
Jako zadanie pozostawiamy policzenie błędów dla innych funkcji, węzłów i typów splajnów kubicznych interpolacyjnych.