Interpolacja trygonometryczna

Załóżmy, że znamy wartości pewnej ciągłej i okresowej funkcji f(x) o okresie 2π w 2n+1 węzłach. Jako bazę interpolacji przyjmujemy zbiór funkcji trygonometrycznych

Otrzymujemy zatem wielomian interpolacyjny w postaci

zawierający 2n+1 nieznanych parametrów.

Ze względu na uproszczenie obliczeń najistotniejszy jest przypadek interpolacji funkcji określonej na zbiorze równoodległych węzłów x_i∈[0,2π] dobranych według następującej zależności

, gdzie i = 0, 1, ..., 2n

Czyli

Warunek interpolacji prowadzi do układu równań liniowych w postaci

Współczynniki pierwszego wiersza macierzy X wynikają z wartości funkcji sin(hx) i cos (hx) dla x₀ = 0. Przedstawiony układ równań rozwiązuje się natychmiastowo, ponieważ macierz X^-1 można wyznaczyć z zależności

Przykład: Daną funkcję f(x) = 7x – x² przybliżyć wielomianem trygonometrycznym przyjmując n = 3. Współrzędne węzłów interpolacji obliczamy ze wzoru:

x	0	0,898	1,795	2,693	3,590	4,488	5,386
y	0	5,478	9,344	11,598	12,242	11,274	8,695

Baza interpolacji – zbiór funkcji:

Interpolacja wielomianowa

W praktyce często używa się bazy złożonej z jednomianów

Baza dla funkcji ciągłych na odcinku skończonym [x_o, x_n] jest bazą zamkniętą, tzn., że każda funkcja tej klasy może być przedstawiona w postaci szeregu złożonego z funkcji bazowych. Wielomian interpolacyjny ma w tym przypadku postać:

dodatkowo musi być spełniony warunek:

Powyższy układ równań ma jedyne rozwiązanie względem a₁, jeżeli wartości x₀, x₁, ..., x_n są od siebie różne

Macierz X^-1 dla bazy wielomianowej nazywana jest macierzą Lagrange’a

Zauważyć należy, że każdy zbiór węzłów równoodległych x_i+1 – x_i = h = const. można sprowadzić do zbioru podstawowego podstawiając , wówczas , a macierz Φ przyjmuje postać

Interpolacja Lagrange’a

Przedstawiony powyżej sposób podejścia do interpolacji nie jest zbyt efektywny, ponieważ macierz X jest macierzą pełną i nie zawsze dobrze uwarunkowaną, co oznacza, że numeryczna procedura jej odwracania może być obarczona dużym błędem.

W interpolacji wielomianowej Lagrange’a dla n+1 węzłów interpolacji

przyjmuje się funkcje bazowe w postaci

Funkcje te są wielomianami stopnia n zbudowanymi w ten sposób, że w funkcji bazowej φ₁ brakuje czynnika (x-x_i). Zatem wielomian interpolacji wyraża się następującym wzorem:

współczynniki a₀ ... a_n tego wielomianu wyznaczamy z równania:

X · A = Y, przy czym macierz X ma postać:

Macierz posiada tylko główną przekątną niezerową w związku z tym układ równań X · A = Y rozwiązuje się natychmiastowo

Można więc wielomian interpolacyjny Lagrange’a zapisać w postaci ułamka:

lub krócej

, j = 0, 1, ..., n

Aproksymacja średniokwadratowa

Niech dana będzie funkcja y=f(x), która w pewnym zbiorze X punktów x₀, x₁, ..., x_n przyjmuje wartości y₀, y₁, ..., y_n. Wartości te znane tylko w przybliżeniu z pewnym błędem (np. jako wyniki pomiarów). Poszukujemy takiej funkcji Q(x) przybliżającej daną funkcję f(x), która umożliwi wygładzenie funkcji f(x), czyli pozwoli z zakłóconych błędami danymi wartości funkcji przybliżonej otrzymać gładką funkcję przybliżającą.

Niech ϕ_j(x), j=0, 1,...,n będzie układem funkcji bazowych. Poszukujemy wielomianu uogólnionego Q(x) będącego najlepszym przybliżeniem średniokwadratowym funkcji f(x) na zbiorze X=(x_j). Funkcja przybliżająca ma postać

Przy czym współczynniki a_i są tak określone, aby wyrażenie

dla i=0, 1, ..., n

było minimalne. Funkcja w(x) jest z góry ustaloną funkcją wagową.

Aby wyznaczyć współczynniki a_i oznaczamy odchylenie

gdzie R_j jest odchyleniem w punkcie x_j.

Następnie obliczamy pochodne cząstkowe funkcji H względem a_i. Z warunku

, gdzie k = 0, 1, ..., n

otrzymujemy układ m+1 równań o niewiadomych a_i zwany układem normalnym

, gdzie k = 0, 1, ..., n

Jeśli wyznacznik tego układu jest różny od zera to rozwiązaniem układu jest minimum funkcji H. W zapisie macierzowym układ przyjmuje postać

gdzie

Aproksymacja wielomianowa

Jeżeli za funkcje bazowe przyjmiemy ciąg jednomianów (xⁱ), i = 0, 1, ..., n to układ normalny przyjmuje postać

, k = 0, 1, ..., n

co po zmianie kolejności sumowania daje

Aproksymacja za pomocą wielomianów ortogonalnych ze wzrostem stopnia wielomianu, obliczenia stają się coraz bardziej pracochłonne a ich wyniki niepewne. Problem ten można usunąć stosując do aproksymacji wielomiany ortogonalne.

Funkcje f(x) i g(x) nazywa się ortogonalnymi na dyskretnym zbiorze punktów x₀, x₁, ..., x_n jeśli

przy czym

Analogicznie ciąg funkcyjny

nazywamy ortogonalnym na zbiorze punktów x₀, x₁, ..., x_n jeśli

dla j≠k

Zastosowanie tej metody powoduje, że znika jedna z trudności obliczeniowych przy aproksymacji wielomianowej, mianowicie złe uwarunkowanie macierzy układu normalnego. Przy aproksymacji wielomianami ortogonalnymi macierz układu normalnego jest macierzą diagonalną, a jej elementy położone na głównej przekątnej dane są wzorem

Załóżmy, że znamy n+1 równoodległych punktów x_i (x_i = x₀+ih, i = 0, 1, ..., n). Za pomocą przekształcenia liniowego

przeprowadzimy te punkty w kolejne liczby całkowite od 0 do n poszukujemy ciągu wielomianów

(dolny indeks oznacza stopień i m ≤ n) spełniających warunek ortogonalności

dla j≠k

przy czym

gdzie

Często używa się unormowanego ciągu wielomianów spełniających warunek

, gdzie k = 0, 1, ..., m

Co po przekształceniu daje nam wzór na wielomiany Grama, zwane też wielomianami Czebyszewa stopni k = 0, 1, ..., m w postaci

, gdzie k = 0, 1, ..., m

Wzór aproksymacyjny oparty na wielomianach Grama ma postać:

gdzie,

oraz

Aproksymacja trygonometryczna

Często spotykamy się z przypadkiem, gdy funkcja f(x) jest okresowa. Taką funkcję wygodniej jest aproksymować, wielomianem trygonometrycznym o bazie:

1, sin x, cos x, sin 2x, cos 2x, …, sin kx, cos kx

Jeżeli f(x) jest funkcją dyskretną określoną w dyskretnym zbiorze równoodległych punktów i ich liczba jest parzysta i wynosi 2L (dla nieparzystej liczby punktów rozumowanie jest analogiczne), niech:

dla i = 0, 1, …, 2L-1

Baza jest ortogonalna nie tylko na przedziale <0, 2π>, ale też na zbiorze punktów x_i, przy czym warunki ortogonalności mają postać:

dla m, k dowolnych, przy czym m i k zmieniają się od 0 do L

Przybliżeniem funkcji f(x) na zbiorze punktów x_i jest wielomian trygonometryczny:

, n<L

Współczynniki a_j i b_j wielomianu wyznaczamy tak, aby suma kwadratów różnic

była minimalna.

Korzystając z warunku ortogonalności otrzymujemy rozwiązanie układu normalnego

w postaci:

dla j=1, 2, …, n

Kwadratury Newtona-Cotesa – wzory zamknięte

Poniżej znajdują się wzory zamknięte dla n = (1..3):

wzór trapezów

wzór parabol (Simpsona)

wzór Bessela

Interpretacja graficzna

wzór parabol (Simpsona)

wzór Bessela

Kwadratury Newtona - Cotesa - wzory otwarte

Wzory kwadratur nie opartych na węzłach będących końcami przedziału całkowania nazywamy wzorami otwartymi Newtona - Cotesa

Jeżeli dla skończonego [a ; b] przedziału wybierzemy zbiór punktów węzłowych {x₀, …, x_n} takich, że :

gdzie:

dla i=(0, …, n)

oraz

W oparciu o Twierdzenie 1 można również wyznaczyć wzory kwadratur otwartych. Wówczas wzory kwadratur dla n parzystych mają postać

oraz dla n nieparzystych mają postać:

wzór prostokątów

Metody Rungego-Kutty

Powszechnie na całym świecie stosuje się metody Rungego–Kutty czwartego rzędu. Polegają one na rozwiązaniu zagadnienia:

gdzie t∈ [a,b] oraz y(a)=α (25)

i przedstawieniu różnicy funkcji y(t) w punktach t_i+1 oraz t_i w postaci:

w_i+1 - w_i = lub inaczej w_i+1 – w_i = h0 (t_i,w_i,h) (26)

gdzie m oznacza rząd metody, c_j są stałymi, a

gdzie α_j, β_j, γ_j, _δlj – stałe

h - krok całkowania

Określenie stałych c_j, α_j, β_j otrzymujemy poprzez rozwinięcie funkcji f(t,y) w szereg Taylora w otoczeniu punktu t_i

Metody R-K – metoda 2 rzędu

Metoda wyprowadzona przez rozwinięcie f(t,y) w szereg Taylora 2 rzędu pozwala określić stałe c₁, α₁,β₁, c₂, α₂, β₂ :

Metoda punktu środkowego:

C1 = 0 C2 = 1 α₁ = 0 α₂ = h/2 β1 = 0 β2 = ½ (28)

wówczas możemy zapisać:

k₁=hf (t_i,w_i) k₂=hf (t_i + ½ h, w_i + ½ k₁)

ostatecznie:

w_i+1=w_i+k₂

lub

w_i+1=w_i + hf(t_i + h/2, w_i + h/2 * f(t_i, w_i))

Interpretacja graficzna:

f₁=f(t_i,w_i) f₂=f(t_i + h/2, w_i + h/2 * f₁)

Stosując metody Rungego–Kutty 2 rzędu , rozwiązać następujące równanie różniczkowe:

y’=y-t² + 1 0≤ t ≤2 y(0)=0,5

dla N = 10 wtedy h = 0,2; t_i = i⋅h ; w₀ = 0,5 oraz

metoda punktu środkowego w_i+1 = 1,22 w_i – 0,0088 i² + 0,218

Metoda zmodyfikowana Eulera

w_i+1 = 1,22 w_i - 0,0088 i² + 0,216

Metoda Heana dla i = 0, 1, ..., 9

w_i+1 = 1,22 w_i - 0,0088 i² + 0,2173

Rozwiązanie dokładne ma postać:

y(t) = (t²+1) - 0,5 e^t

Metody R-K rzędu 4

Metoda wyprowadzona przez rozwinięcie f(t,y) w szereg Taylora 4 rzędu, pozwala określić stałe we wzorze (wzór ogólny R-K). Poniżej przedstawiono najczęściej stosowaną metodę 4 rzędu

k₁ = hf (t_i, w_i)

k₂ = hf (t_i + ½ *h, w_i + ½ *k₁)

k₃ = hf (t_i + ½ *h, w_i + ½ *k₂) (37)

k₄ = hf (t_i+h, w_i+k₃)

ostatecznie:

w_i+1=w_i + 1/6 (k₁ + 2k₂ + 2k₃ + k₄) (38)

Interpretacja graficzna:

f₁ = f (t_i, w_i)

f₂ = f (t_i + ½ * h, w_i+1 + h f₁)

f₃ = f (t_i + ½ * h, w_i + ½ * h f₂)

f₄ = f (t_i + h, w_i + h f₃)

linia 3-4 nie jest łamana – jest prosta !!! (na 99%) :)

Metoda R-K jest najczęściej stosowaną metodą do rozwiązywania układów równań różniczkowych

Tworzymy macierz X:

Elementy macierzy X mnożymy przez 2/7, otrzymaną macierz transponujemy i obliczamy współczynniki wzoru interpolacyjnego

Metoda Eulera

Rozpatrzmy zadanie znalezienia funkcji y=y(t), które dla spełnia równanie różniczkowe zwyczajne pierwszego rzędu:

(7)

oraz warunek początkowy

y(a)=y₀ (8),

gdzie f jest daną funkcją dwu zmiennych. Warunek wystarczający istnienia i jednoznaczności rozwiązania zadania (7) oraz (8) jest warunek ciągłości funkcji f(t,y)

Metoda ta rzadko jest używana w praktyce, natomiast ze względu na jej prostotę, zostanie przedstawiona w celu zobrazowania i zrozumienia technik stosowanych w bardziej zaawansowanych metodach, pomijając przy tym złożone działania matematyczne. Przedmiotem rozważań będzie poniższe równanie.

, gdzie oraz y(a)=a (10)

Na rozwiązanie powyższego zagadnienia będziemy obliczać przybliżone wartości funkcji yi=y(ti), gdzie ti=a+ih, h=(b-a)/N, dla którego i=(0,1, ..., N), gdzie ti nazywane są punktami węzłowymi, natomiast h odległością między nimi.

Rozwiniemy y(t) w szereg Taylora w celu wyprowadzenia metody Eulera. Zakładając, że y(t) jest jedynym rozwiązaniem (10) oraz posiada drugą pochodną, która jest ciągła w przedziale [a, b] wówczas dla każdego i=1, 2, ..., N.

Wykorzystując powyższe założenia możemy zapisać:

(11), gdzie

Oznaczamy h=t_i+1-t_i, wówczas otrzymujemy:

(12)

Użyjemy podstawienia y’(t)=f(t, y), wówczas otrzymujemy:

(13)

Zapisując ω_i≈y(t_i) oraz pomijając błąd przybliżenia otrzymujemy ω_i+1=ω_i+hf(t_i, ω_i) (14)

Powyższy wzór nazywamy metodą Eulera – wzór ten nazywany jest inaczej równaniem różniczkowym, gdyż można zapisać:

(15)

Aby wyznaczyć wartość szukanej funkcji y(x) w następnym kroku h, wykorzystujemy poprzednią wartość funkcji oraz wielkości zmian funkcji – dzięki pochodnej. Natomiast uwzględniając błąd przybliżenia wzór (13) przyjmuje postać:

, gdzie (16)

Lokalny błąd dyskretyzacji τ_i+1(h)

(17)

dodatkowo możemy określić krok h, dla którego błąd lokalny jest mniejszy od zadanej dokładności δ

, gdzie

Globalny błąd dyskretyzacji g(x)

g(t)= ω(t)-y(t) (19)

Dla wybranego punktu t_i możemy zapisać:

g_i=ω_i-y(t_i) (20), wówczas

, gdzie L- liczba Lipschnitz`a