Metody numeryczne: stycznia 2014

Kwadratura Gaussa

Kwadraturami Gaussa nazywamy metody całkowania numerycznego polegające na takim wyborze wag $w_1, w_2, \ldots, w_n$ i węzłów interpolacji $t_1, t_2, \ldots, t_n \in [a,b]$ aby wyrażenie

$\sum_{i=1}^n w_i f(t_i)$

najlepiej przybliżało całkę

$I(f) = \int\limits_a^b w(x) f(x) dx$

gdzie $f$ jest dowolną funkcją określoną na odcinku $[a,b]$ , a $w$ jest tzw. funkcją wagową spełniającą warunki

$w(x)\ge 0$ ,
$\forall_{k\in \mathbb{N}} \int\limits_a^b x^k w(x) dx$ jest skończona,
Jeżeli $p$ jest wielomianem takim, że $\forall_{x\in [a,b]}\;p(x)\ge 0$ , to jeśli $\int\limits_a^b w(x) p(x) dx=0$ , mamy wtedy $p\equiv0$ .

Określmy iloczyn skalarny z wagą

$\langle f,g \rangle_w = \int\limits_a^b w(x) f(x) g(x) dx.$

Powiemy, że dwa wielomiany są ortogonalne względem tego iloczynu skalarnego jeśli $\langle f,g \rangle_w =0$ .

Wszystkie kwadratury Gaussa wywodzą się z twierdzenia udowodnionego przez niego:

a) Jeżeli $t_1, t_2, \ldots, t_n \in [a,b]$ są pierwiastkami n-tego wielomianu ortogonalnego $p_n(x)$ oraz $w_1, w_2, \ldots, w_n$ są rozwiązaniami układu równań:

$\left\{ \begin{matrix} p_0(t_1)w_1 + & \ldots & + p_0(t_n)w_n= & <p_0,p_0>_w\\ p_1(t_1)w_1 + & \ldots & + p_1(t_n)w_n= & 0\\ \vdots & & \vdots & \vdots \\ p_{n-1}(t_1)w_1 + & \ldots & + p_{n-1}(t_n)w_n= & 0\\ \end{matrix}\right.$

to dla każdego wielomianu $p$ stopnia nie większego niż 2n-1 zachodzi

$\int\limits_a^b w(x) p(x) dx = \sum_{i=1}^n w_i p(t_i).\qquad (*)$

Ponadto $w_i > 0$ .

b) Jeżeli dla pewnego ciągu węzłów $x_1, x_2, \ldots, x_n \in [a,b]$ oraz ciągu wag $v_1, v_2, \ldots, v_n$ dla dowolnego wielomianu $p$ stopnia nie większego niż 2n-1 zachodzi warunek (*), to $x_i=t_i$ oraz $v_i=w_i$ z dokładnością do kolejności.

c) Dla dowolnego ciągu węzłów $x_1, x_2, \ldots, x_n \in [a,b]$ oraz ciągu wag $v_1, v_2, \ldots, v_n$ nie istnieje wielomian stopnia 2n, dla którego nie zachodzi warunek (*).

Wzór trapezów

Wzór trapezów – jeden z wielu wzorów służących do przybliżonego obliczania całek oznaczonych w sensie Riemanna. Idea wzoru opiera się na geometrycznej interpretacji całki oznaczonej z funkcji nieujemnej jako pola pod wykresem funkcji.

Jeżeli przedział całkowania [a, b] podzielony zostanie punktami x₁, x₂, ..., x_n-1 na n równych części o długościach (b-a)/n, i w figurę ograniczoną na prostymi x = a, x = b, osią odciętych oraz wykresem funkcji y = f(x) wpiszemy trapezy jak pokazano na rysunku poniżej,

to pola kolejnych trapezów wynoszą:

$\frac{b-a}{n}\cdot\frac{f(x_0)+f(x_1)}{2},\ \frac{b-a}{n}\cdot\frac{f(x_1)+f(x_2)}{2},\ \dots, \frac{b-a}{n}\cdot\frac{f(x_{n-1})+f(x_n)}{2}$

gdzie dla jednolitości oznaczono a = x₀ i b = x_n.

Suma pól trapezów jest w przybliżeniu równa polu całego obszaru, czyli:

$\int\limits_a^b f(x)\,dx\approx\frac{b-a}{{2}{n}}\left(f(x_0)+f(x_n)+2f(x_1)+\dots+2f(x_{n-1})\right)$ .

Ten właśnie wzór nazywany jest wzorem trapezów.

W przypadku funkcji ciągłej na przedziale [a, b], wzór trapezów pozwala obliczać jej całkę oznaczoną na tym przedziale z dowolną dokładnością, wystarczy w tym celu wziąć za n odpowiednio dużą liczbę. Błąd przybliżenia daje się oszacować w przypadku funkcji, która ma na przedziale [a, b] ciągłą drugą pochodną:

$|R_n|\le\left|\frac{(b-a)^3}{12n^2}K\right|$

gdzie K oznacza największą wartość funkcji |f ′′(x)| w przedziale [a, b].

Obecnie wzór trapezów ma znaczenie wyłącznie historyczne – dostępne programy do całkowania numerycznego stosują o wiele dokładniejsze metody i pozwalają uniknąć czasochłonnych rachunków.

Fourier

Szybka transformacja Fouriera (ang. FFT od Fast Fourier Transform) to algorytm liczenia dyskretnej transformaty Fouriera oraz transformaty do niej odwrotnej.

Czasem używana jest też forma szybka transformata Fouriera w odniesieniu do tej metody. Ściśle jednak transformacja jest przekształceniem, a transformata wynikiem tego przekształcenia.

Niech x₀, ...., x_N-1 będą liczbami zespolonymi, wtedy dyskretna transformata Fouriera jest określona wzorem

$X_k = \sum_{n=0}^{N-1} x_n e^{-{2\pi i \over N} nk } \qquad k = 0,\dots,N-1.$

Obliczanie tych sum za pomocą powyższego wzoru zajęłoby O(N²) operacji.

Algorytmy (jak algorytm Cooleya-Tukeya) obliczające szybką transformację Fouriera bazują na metodzie dziel i zwyciężaj rekurencyjnie dzieląc transformatę wielkości N = N₁N₂ na transformaty wielkości N₁ i N₂ z wykorzystaniem O(N) operacji mnożenia.

Najpopularniejszą wersją FFT jest FFT o podstawie 2. Jest to bardzo efektywna operacja, jednak wektor próbek wejściowych (spróbkowany sygnał) musi mieć długość $N=2^k$ , gdzie $k$ to pewna liczba naturalna. Wynik otrzymuje się na drodze schematycznych przekształceń, opartych o tak zwane struktury motylkowe.

Złożoność obliczeniowa Szybkiej transformacji Fouriera wynosi $O(N \log_2 N)$ , zamiast $O(N^2)$ algorytmu wynikającego wprost ze wzoru określającego dyskretną transformatę Fouriera. Dzięki istnieniu takiego algorytmu praktycznie możliwe stało się cyfrowe przetwarzanie sygnałów (DSP), a także zastosowanie dyskretnych transformat kosinusowych (DCT) (JPEG,MP3 itd.) do kompresji.

Monte Carlo

Metoda Monte Carlo (MC) jest stosowana do modelowania matematycznego procesów zbyt złożonych (obliczania całek, łańcuchów procesów statystycznych), aby można było przewidzieć ich wyniki za pomocą podejścia analitycznego. Istotną rolę w metodzie MC odgrywa losowanie (wybór przypadkowy) wielkości charakteryzujących proces, przy czym losowanie dokonywane jest zgodnie z rozkładem, który musi być znany.

Typowym przykładem może być modelowanie wyniku zderzenia cząstki o wysokiej energii z jądrem złożonym, gdzie każdy akt zderzenia elementarnego (z pojedynczym nukleonem jądra) modelowany jest oddzielnie poprzez losowanie liczby, rodzaju, kąta emisji, energii itp. cząstek wtórnych emitowanych w wyniku takiego zderzenia. Następnym etapem jest modelowanie losu każdej z cząstek wtórnych (w wyniku kolejnego losowania prawdopodobieństwa oddziaływania lub wyjścia z jądra). Kontynuując taką procedurę można otrzymać pełny opis "sztucznie generowanego" procesu złożonego. Po zebraniu dostatecznie dużej liczby takich informacji można zestawić ich charakterystyki z obserwowanymi wynikami doświadczalnymi, potwierdzając lub negując słuszność poczynionych w całej procedurze założeń.

Metoda została opracowana i pierwszy raz zastosowana przez Stanisława Ulama.

Metodą Monte Carlo można obliczyć pole figury zdefiniowanej nierównością:

$x^2 +y^2 \le R^2$

czyli koła o promieniu R i środku w punkcie (0,0).

Losuje się n punktów z opisanego na tym kole kwadratu - dla koła o R = 1 współrzędne wierzchołków (-1,-1), (-1,1), (1,1), (1,-1).
Po wylosowaniu każdego z tych punktów trzeba sprawdzić, czy jego współrzędne spełniają powyższą nierówność (tj. czy punkt należy do koła).

Wynikiem losowania jest informacja, że z n wszystkich prób k było trafionych, zatem pole koła wynosi

$P_k=P\ \frac {k}{n}$

gdzie P jest polem kwadratu opisanego na kole (Dla R = 1 : P = 4).

Wzory Newtona-Cotesa

W analizie numerycznej wzory Newtona-Cotesa są zbiorem metod numerycznych całkowania, zwanego również kwadraturą. Nazwa pochodzi od Isaaca Newtona i Rogera Cotesa.

Przyjmujemy, że wartość funkcji f jest znana w równo oddalonych punktach x_i, dla i = 0, ..., n. Dla punktów oddalonych od siebie o inne odległości ma zastosowanie inna klasa wzorów,kwadratura gaussowska.

Jeżeli $a=x_0<x_1<x_2,\cdots <x_{n-1}<x_n=b$ są równo odległymi węzłami interpolacji funkcji f(x) (tj. $x_i$ są elementami dziedziny f, dla których znana jest wartość $f(x_{i })=y_i$ ), to całkę:

$\int\limits_a^b f(x) dx$

można aproksymować całką:

$\int\limits_a^b L_{n }(x) dx$

gdzie $L_{n }(x) dx$ jest wielomianem interpolacyjnym Lagrange'a stopnia co najwyżej n, przybliżającym funkcję f(x) w węzłach interpolacji, tj.:

$L_{n }(x_{0 })=y(x_{0 }),L_{n }(x_{1 })=y(x_{1 }),\cdots ,L_{n }(x_{n })=y(x_{n })$

Niech $h_n = \frac {b-a}{n}$ oznacza długość kroku dzielącą dwa węzły interpolacji.

Wprowadzając zmienną t taką, że x=a+th można zapisać:

$L_{i }(x)=L_{n }(a+th)=\prod_{j = 0 \and j\ne i}^n \frac {t-j}{i-j}=g(t)$

Wtedy:

$\int\limits_a^b L_{n }(x) dx=\int\limits_a^b \sum_{i = 0}^n f(x_i)\cdot L_{i }(x) dx=$

$\sum_{i = 0}^n f(x_i)\cdot \int\limits_a^b L_{i }(x) dx$

$x=a+t\cdot h$ ,

$f(x_i)=f(a+i\cdot h)$

$x_i=a+i\cdot h$

dla $a=x_0$ t=0

dla $x_1$ t=1

$\cdots$

dla $b=x_n$ t=n

$dx=(x)'=1$

$dt=(a+t\cdot h)dt=(a+t\cdot h)'=h=dt$

Zmieniając zmienną, oraz granice całkowania otrzymuje się:

$\int\limits_a^b L_{i }(x)dx = h\cdot \int\limits_0^n g(t) dt$

Ostatecznie, wzór Newtona-Cotesa dla n+1 równo odległych węzłów przyjmuje postać:

$\int\limits_a^b f(x) dx = \int\limits_a^b L_{n }(x) dx = \sum_{i=0}^n f(x_{i }) \cdot \int\limits_a^b L_{i }(x=a+t\cdot h) dx = \sum_{i=0}^n f(x_{i }) \cdot h\cdot \int\limits_0^n \prod_{j=0 \and j\ne i}^n \frac {t-j}{i-j} dt$

Przyjmując za $A_i = h\cdot \int\limits_0^n \prod_{j=0 \and j\ne i}^n \frac {t-j}{i-j} dt$ (nazywane współczynnikami kwadratury Newtona-Cotesa), otrzymuje się: