MN14

From Studia Informatyczne


Spis treści

Całkowanie

<<< Powrót do strony głównej przedmiotu Metody numeryczne

Zajmiemy się teraz zadaniem całkowania numerycznego. Polega ono na obliczeniu (a raczej przybliżeniu) całki oznaczonej

\displaystyle S(f)\,=\,\int_a^b f(x)\,dx,

gdzie \displaystyle -\infty<a<b<+\infty, a \displaystyle f należy do pewnej klasy \displaystyle F funkcji rzeczywistych określonych i całkowalnych w sensie Riemanna na całym przedziale \displaystyle [a,b].

Każdy, kto przeszedł przez kurs całkowania wie, że obliczanie całek rozumiane jako znalezienie elementarnego wzoru na funkcję pierwotną może być trudne, bardzo trudne, a nawet niewykonalne. Tymczasem zadanie przybliżonego wyznaczenia wartości całki daje się w dużej mierze zautomatyzować z całkiem dobrym skutkiem.

Obliczanie całek jest wymagane w bardzo wielu zadaniach inżynierskich i naukowych. Całki z funkcji (bardzo) wielu zmiennych (które na swój sposób są szczególnie trudne do obliczenia) znajdują ważne zastosowania w bankowości i finansach.

Będziemy zakładać, że mamy możliwość obliczania wartości funkcji \displaystyle f, a w niektórych przypadkach również jej pochodnych, o ile istnieją. Dokładna całka \displaystyle S(f) będzie więc w ogólności przybliżana wartością \displaystyle A(f), która zależy tylko od wartości \displaystyle f i ewentualnie jej pochodnych w skończonej liczbie punktów.

Kwadratury

Kwadraturami nazywamy funkcjonały liniowe \displaystyle Q:F\to R postaci

\displaystyle Q(f)\,=\,\sum_{i=0}^n a_i f(x_i),

albo ogólniej

\displaystyle    Q(f)\,=\,\sum_{i=0}^k\sum_{j=0}^{n_i-1}                   a_{i,j}f^{(j)}(x_i),

gdzie \displaystyle x_i są punktami z \displaystyle [a,b], a \displaystyle a_i (albo \displaystyle a_{i,j}) są pewnymi współczynnikami rzeczywistymi. Zauważmy, że obliczenia kwadratur są dopuszczalne w naszym modelu obliczeniowym, mogą więc służyć jako sposób przybliżania całki.

Jeden z możliwych sposobów konstrukcji kwadratur jest następujący. Najpierw wybieramy węzły \displaystyle x_j (pojedyncze lub wielokrotne), budujemy wielomian interpolacyjny odpowiadający tym węzłom, a następnie całkujemy go. Ponieważ postać wielomianu interpolacyjnego zależy tylko od danej informacji o \displaystyle f, otrzymana w ten sposób wartość też będzie zależeć tylko od tej informacji, a w konsekwencji funkcjonał wynikowy będzie takiej postaci, jak wyżej. Są to tzw. kwadratury interpolacyjne.

Definicja

Kwadraturę \displaystyle Q^I opartą na węzłach o łącznej krotności \displaystyle n+1 nazywamy interpolacyjną, jeśli

\displaystyle Q^{I}(f)\,=\,\int_a^b w_f(x)\,dx,

gdzie \displaystyle w_f jest wielomianem interpolacyjnym funkcji \displaystyle f stopnia co najwyżej \displaystyle n, opartym na tych węzłach.

Współczynniki kwadratur interpolacyjnych można łatwo wyliczyć. Rozpatrzmy dla uproszczenia przypadek, gdy węzły są jednokrotne. Zapisując wielomian interpolacyjny w postaci jego rozwinięcia w bazie kanonicznej Lagrange'a \displaystyle l_i, otrzymujemy

\displaystyle Q^{I}(f) \,=\, \int_a^b \sum_{i=0}^n f(x_i)l_i(x)\,dx             \,=\, \sum_{i=0}^n f(x_i)\int_a^b l_i(x)\,dx,

a stąd i z postaci \displaystyle l_i,

\displaystyle a_i\,=\,\int_a^b \frac    {(x-x_0)\cdots(x-x_{i-1})(x-x_{i+1})\cdots(x-x_n)}    {(x_i-x_0)\cdots(x_i-x_{i-1})(x_i-x_{i+1})\cdots(x_i-x_n)}     \,dx,

\displaystyle 0\le i\le n.

Podamy teraz kilka przykładów.

Kwadratura prostokątów jest oparta na jednym węźle \displaystyle x_0=(a+b)/2,

\displaystyle Q^{I}_0(f)\,=\,(b-a)f\Big(\frac{a+b}2\Big).
Kwadratura prostokątów
Enlarge
Kwadratura prostokątów

Kwadratura trapezów jest oparta na jednokrotnych węzłach \displaystyle x_0=a, \displaystyle x_1=b i jest równa polu odpowiedniego trapezu,

\displaystyle Q^{I}_1(f)\,=\,T(f)\,=\,\frac{b-a}2 \Big(f(a)+f(b)\Big).
Kwadratura trapezów
Enlarge
Kwadratura trapezów

Kwadratura parabol (Simpsona) jest oparta na jednokrotnych węzłach \displaystyle x_0=a, \displaystyle x_1=b, \displaystyle x_2=(a+b)/2, i jest równa polu pod parabolą interpolującą \displaystyle f w tych węzłach,

\displaystyle Q^I_2(f)\,=\,P(f)\,=\,\frac{b-a}6     \Big( f(a)+4f\Big(\frac{a+b}2\Big)+f(b) \Big).

Zauważmy, że kwadratury trapezów i parabol są oparte na węzłach jednokrotnych i równoodległych, przy czym \displaystyle x_0=a i \displaystyle x_n=b. Ogólnie, kwadratury interpolacyjne oparte na węzłach równoodległych \displaystyle x_i=a+(b-a)i/n, \displaystyle 0\le i\le n, nazywamy kwadraturami Newtona--Cotesa.

Błąd kwadratur interpolacyjnych

Zajmiemy się teraz błędem kwadratur interpolacyjnych. Przypomnijmy, że \displaystyle F^r_M([a,b]) oznacza klasę funkcji \displaystyle (r+1) razy różniczkowalnych w sposób ciągły i takich, że \displaystyle |f^{(r+1)}(x)|\le M, \displaystyle \forall x.

Twierdzenie O błędzie kwadratur interpolacyjnych

Niech \displaystyle Q^I będzie kwadraturą interpolacyjną opartą na (jednokrotnych lub wielokrotnych) węzłach \displaystyle x_i, \displaystyle 0\le i\le n. Jeśli \displaystyle f\in F^n_M([a,b]), to

\displaystyle    |S(f)\,-\,Q^I(f)|\,\le\,\frac M{(n+1)!}(b-a)^{n+2}.

W klasie \displaystyle F^n_M([a,b]) maksymalny błąd kwadratury \displaystyle Q^I wynosi

\displaystyle \sup_{f\in F_M^n([a,b])} |S(f)-Q^I(f)|\,=\,\frac M{(n+1)!}      \int_a^b |(x-x_0)(x-x_1)\cdots(x-x_n)|\,dx.

Dowód

Korzystając ze znanego nam już wzoru na błąd interpolacji wielomianowej, mamy

\displaystyle S(f)\,-\,Q^{I}(f)\,=\,\int_a^b   (x-x_0)(x-x_1)\cdots(x-x_n)f(x_0,x_1,\ldots,x_n,x)\,dx.

Stąd, jeśli \displaystyle f\in F^n_M([a,b]), to

\displaystyle |S(f)\,-\,Q^I(f)|\,\le\,\int_a^b (b-a)^{n+1}\frac M{(n+1)!}\,dx      \,=\,(b-a)^{n+2}\frac M{(n+1)!}.

Ograniczenie górne w dokładnej formule na błąd w klasie \displaystyle F^n_M([a,b]) wynika bezpośrednio. Aby pokazać ograniczenie dolne zauważmy, że dla funkcji \displaystyle g takiej, że \displaystyle g^{(n+1)} przyjmuje na przedziałach \displaystyle (a,x_0), \displaystyle (x_0,x_1), \displaystyle \ldots, \displaystyle (x_n,b) naprzemiennie wartości \displaystyle M i \displaystyle -M mamy

\displaystyle |S(g)-Q^I(g)|\,=\,\frac M{(n+1)!}      \int_a^b |(x-x_0)(x-x_1)\cdots(x-x_n)|\,dx.

Co prawda, \displaystyle g nie jest w \displaystyle F^n_M([a,b]), ale może być dla dowolnego \displaystyle \epsilon>0 przybliżana funkcjami \displaystyle f_\epsilon\in F^n_M([a,b]) w ten sposób, że całka

\displaystyle \int_a^b |(x-x_0)\cdots(x-x_n)(f-g)^{(n+1)}(x)|\,dx     \,\le\,\epsilon\,(n+1)!.

Zapisując \displaystyle f_\epsilon=g+(f_\epsilon-g) mamy

\displaystyle \aligned |S(f_\epsilon)\,-\,Q^I(f_\epsilon)| &\le & |S(g)\,-\,Q^I(g)|\,+\,       |S(f_\epsilon-g)-Q^I(f_\epsilon-g)| \\     &\le & \frac M{(n+1)!}\int_a^b |(x-x_0)\cdots(x-x_n)|      \,dx\,+\,\epsilon,  \endaligned

co wobec dowolności \displaystyle \epsilon daje dowód twierdzenia.

image:End_of_proof.gif

W szczególnych przypadkach kwadratur trapezów \displaystyle T i parabol \displaystyle P możemy otrzymać innego rodzaju formuły na błąd.

Twierdzenie O postaci błędu kwadratury trapezów i Simpsona

Jeśli \displaystyle f\in C^{(2)}([a,b]), to dla kwadratury trapezów mamy

\displaystyle S(f)\,-\,T(f)\,=\,-\frac{(b-a)^3}{12}f^{(2)}(\xi_1).

Jeśli \displaystyle f\in C^{(4)}([a,b]), to dla kwadratury parabol mamy

\displaystyle S(f)\,-\,P(f)\,=\,-\frac{(b-a)^5}{2280}f^{(4)}(\xi_2).

(\displaystyle \xi_1,\xi_2\in [a,b]).

Dowód

Najpierw udowodnimy część dotyczącą kwadratury trapezów. Ze wzoru na błąd kwadratury,

\displaystyle S(f)\,-\,T(f)\,=\,\int_a^b (x-a)(x-b)f(a,b,x)\,dx.

Ponieważ funkcja \displaystyle x\mapsto f(a,b,x) jest ciągła, a wielomian \displaystyle (x-a)(x-b) przyjmuje jedynie wartości nieujemne, można zastosować twierdzenie o wartości średniej dla całki, aby otrzymać

\displaystyle \aligned S(f)\,-\,T(f) &= f(a,b,c)\int_a^b (x-a)(x-b)\,dx \\     &= -\frac{f^{(2)}(\xi_1)}{2!}\frac{(b-a)^3}6, \endaligned

dla pewnych \displaystyle c,\xi_1\in [a,b].

Teraz zajmiemy się kwadraturą parabol. Niech \displaystyle w_{f,2}\in\Pi_2 i \displaystyle w_{f,3}\in\Pi_3 będą wielomianami interpolacyjnymi funkcji \displaystyle f odpowiednio dla węzłów \displaystyle a,b,(a+b)//2 oraz \displaystyle a,b,(a+b)//2,(a+b)//2. Wtedy

\displaystyle w_{f,3}(x)\,=\,w_{f,2}(x)\,+\,      f\Big(a,b,\frac{a+b}2,\frac{a+b}2\Big)        (x-a)\Big(x-\frac{a+b}2\Big)(x-b).

Wobec

\displaystyle \int_a^b (x-a)\Big(\frac{a+b}2\Big)(x-b)\,dx\,=\,0

mamy

\displaystyle P(f) \,=\, \int_a^b w_{f,2}(x)\,dx\,=\,            \int_a^b w_{f,3}(x)\,dx.

Stąd i ze wzoru na błąd interpolacji Hermite'a otrzymujemy

\displaystyle \aligned \lefteqn{ S(f)\,-\,P(f)\;=\;\int_a^b (f-w_{f,3})(x)\,dx } \\     && =\; \int_a^b (x-a)\Big(x-\frac{a+b}2\Big)^2(x-b)         f\Big(a,b,\frac{a+b}2,\frac{a+b}2,x\Big)\,dx.  \endaligned

Ponieważ wielomian \displaystyle (x-a)(x-(a+b)/2)^2(x-b) jest niedodatni na \displaystyle [a,b], możemy znów zastosować twierdzenie o wartości średniej. Mamy

\displaystyle \aligned S(f)\,-\,P(f) &= f\Big(a,b,\frac{a+b}2,\frac{a+b}2,c\Big) \\    && \qquad\qquad\qquad    \int_a^b (x-a)\Big(x-\frac{a+b}2\Big)^2(x-b)\,dx \\    &= -\frac{f^{(4)}(\xi_2)}{4!}\frac{(b-a)^5}{120},  \endaligned

co kończy dowód.

image:End_of_proof.gif

Kwadratury złożone

Chcielibyśmy, aby błąd kwadratur malał do zera, gdy liczba węzłów rośnie do nieskończoności. Można to osiągnąć stosując np. kwadratury złożone. Są to kwadratury, które powstają przez scałkowanie funkcji kawałkami wielomianowej interpolującej \displaystyle f.

Georg Riemann  Zobacz biografię
Enlarge
Georg Riemann
Zobacz biografię

Prostym przykładem kwadratury złożonej jest suma Riemanna,

\displaystyle \bar Q(f)\,=\,\sum_{i=0}^n (t_{i+1}-t_i)f(x_i),

gdzie \displaystyle a=t_0<t_1<\cdots<t_{n+1}=b oraz \displaystyle x_i\in [t_i,t_{i+1}]. Jeśli średnica podziału, \displaystyle \max_{0\le i\le n}(t_i-t_{i-1}), maleje do zera, to \displaystyle \lim_{n\to\infty}\bar Q(f)=S(f).

Będziemy rozpatrywać kwadratury złożone postaci

\displaystyle \bar Q(f)\,=\,\int_a^b \bar w_f(x)\,dx,

gdzie \displaystyle \bar w_f jest kawałkami wielomianem. Dokładniej, dla danego \displaystyle n kładziemy \displaystyle t_i=a+(b-a)i/k, \displaystyle 0\le i\le k, a następnie dla każdego \displaystyle i wybieramy dowolne węzły \displaystyle x_{i,j}\in [t_{i-1},t_i], \displaystyle 0\le j\le r. Wtedy \displaystyle \bar w_f jest na każdym przedziale wielomianem interpolacyjnym funkcji \displaystyle f stopnia co najwyżej \displaystyle r opartym na węzłach \displaystyle x_{i,j}. Kwadratura \displaystyle \bar Q korzysta z węzłów o łącznej krotności \displaystyle n\le k(r+1).

Twierdzenie O błędzie kwadratur złożonych

Błąd kwadratury złożonej \displaystyle \bar Q(f) w klasie \displaystyle F^r_M([a,b]) jest ograniczony przez

\displaystyle \sup_{f\in F^r_M([a,b])} |S(f)-\bar Q(f)|\,\le\,      \frac{(b-a)^{r+2}}{k^{r+1}}      \frac{M}{(r+1)!}\,\le\,C\,\Big(\frac 1n\Big)^{r+1},

gdzie

\displaystyle C\,=\,\frac{M(r+1)^{r+1}(b-a)^{r+2}}{(r+1)!}.

Dowód

Twierdzenie to jest bezpośrednim wnioskiem z twierdzenia o błędzie kwadratur interpolacyjnych. Mamy bowiem

\displaystyle \aligned |S(f)-\bar Q(f)| &\le & \sum_{i=1}^k      \int_{t_{i-1}}^{t_i} |f(x)-\bar w_f(x)|\,dx \\      &\le & \sum_{i=1}^k \Big(\frac{b-a}{k}\Big)^{r+2}               \frac M{(r+1)!} \,=\,      \frac{(b-a)^{r+2}}{k^{r+1}}\frac M{(r+1)!},  \endaligned

co kończy dowód.

image:End_of_proof.gif

W klasie \displaystyle F^r_M([a,b]), błąd kwadratur złożonych jest rzędu \displaystyle n^{-(r+1)}. Można pokazać, że błąd każdej innej metody całkowania korzystającej jedynie z wartości funkcji w \displaystyle n punktach nie może w klasie \displaystyle F^r_M([a,b]) maleć szybciej niż \displaystyle n^{-(r+1)}. Podane kwadratury złożone mają więc optymalny rząd zbieżności.

Zajmiemy się teraz błędem szczególnych kwadratur złożonych, mianowicie złożonych kwadratur trapezów \displaystyle \bar T_k i parabol \displaystyle \bar P_k. Powstają one przez zastosowanie na każdym przedziale \displaystyle [t_{i-1},t_i] odpowiednio kwadratur trapezów \displaystyle T i parabol \displaystyle P.

Złożona kwadratura trapezów
Enlarge
Złożona kwadratura trapezów

Jak łatwo się przekonać,

\displaystyle \bar T_k(f)\,=\,\frac{b-a}{k}\left(\frac{f(a)+f(b)}2     \,+\,\sum_{j=1}^{k-1} f\Big(\frac jk\Big)\right),

oraz

\displaystyle \bar P_k(f)\,=\,\frac{b-a}{3k}\left(\frac{f(a)+f(b)}2     \,+\,\sum_{j=1}^{k-1} f\Big(\frac jk\Big)\,+\,         2\,\sum_{j=1}^k f\Big(\frac{2j-1}{2k}\Big)\right).

Twierdzenie O postaci błędu złożonych kwadratur trapezów i Simpsona

Jeśli \displaystyle f\in C^{(2)}([a,b]), to

\displaystyle S(f)\,-\,\bar T_k(f)\,=\,-\frac{(b-a)^3}{12\,k^2}        f^{(2)}(\xi_1).

Jeśli \displaystyle f\in C^{(4)}([a,b]), to

\displaystyle S(f)\,-\,\bar P_k(f)\,=\,-\frac{(b-a)^5}{2280\,k^4}        f^{(4)}(\xi_2).

Dowód

Dla kwadratury trapezów mamy

\displaystyle \aligned \lefteqn{ S(f)\,-\,\bar T_k(f) \;=\; -\sum_{i=1}^k       \frac{(b-a)^3}{12 k^3}f^{(2)}(\alpha_i) } \\      && =\;-\frac{(b-a)^3}{12 k^2}\frac 1k       \sum_{i=1}^k f^{(2)}(\alpha_i) \,=\,      -\frac{(b-a)^3}{12 k^2} f^{(2)}(\xi_1), \endaligned

a dla kwadratury parabol podobnie

\displaystyle \aligned \lefteqn{ S(f)\,-\,\bar P_k(f) \;=\; -\sum_{i=1}^k       \frac{(b-a)^5}{2280 k^5}f^{(4)}(\beta_i) } \\      && =\; -\frac{(b-a)^5}{2280 k^4}\frac 1k       \sum_{i=1}^k  f^{(4)}(\beta_i) \,=\,      -\frac{(b-a)^5}{2280 k^4} f^{(4)}(\xi_2). \endaligned
image:End_of_proof.gif

Kwadratura parabol ma więc optymalny rząd zbieżności nie tylko w klasie \displaystyle F^2_M([a,b]), ale też w \displaystyle F^3_M([a,b]).

Przyspieszanie zbieżności kwadratur

W praktyce często stosuje się obliczanie kwadratur poprzez zagęszczanie podziału przedziału \displaystyle [a,b]. Na przykład, dla złożonej kwadratury trapezów zachodzi następujący wygodny wzór rekurencyjny:

\displaystyle    \bar T_{2k}\,=\,\frac 12\left(\bar T_k(f)\,+\,       \frac{b-a}k\,\sum_{i=1}^k        f\Big(\frac{2i-1}{2k}\Big)\right).

Pozwala on obliczyć \displaystyle \bar T_{2k}(f) na podstawie \displaystyle \bar T_k(f) poprzez "doliczenie" wartości funkcji w punktach "gęstszej" siatki. W ten sposób możemy obserwować zachowanie się kolejnych przybliżeń \displaystyle \bar T_{2^s}(f) (\displaystyle s\ge 0) całki \displaystyle S(f). Jest to szczególnie istotne wtedy, gdy nie mamy żadnej informacji a priori o \displaystyle \|f''\|_{ C([a,b])}, a przez to nie potrafimy oszacować liczby \displaystyle n węzłów, dla której osiągniemy pożądaną dokładność.

Jeśli funkcja jest więcej niż dwa razy różniczkowalna, to użycie złożonych kwadratur trapezów zdaje się tracić sens. Wtedy istnieją przecież kwadratury, których błąd maleje do zera szybciej niż \displaystyle n^{-2}. Okazuje się jednak, że kwadratury \displaystyle \bar T_k mogą być podstawą dla prostej rekurencyjnej konstrukcji innych kwadratur posiadających już optymalną zbieżność. Konstrukcja ta bazuje na następującym ważnym lemacie.

Leonhard Euler  Zobacz biografię
Enlarge
Leonhard Euler
Zobacz biografię
Colin Maclaurin  Zobacz biografię
Enlarge
Colin Maclaurin
Zobacz biografię

Lemat Formuła Eulera-Maclaurina

Dla funkcji \displaystyle f\in C^{(2m+2)}([a,b]), błąd złożonej kwadratury trapezów \displaystyle \bar T_k wyraża się wzorem

\displaystyle \aligned S(f)\,-\,\bar T_k(f) &= \sum_{i=1}^{m} c_ih^{2i}   \Big(f^{(2i-1)}(b)-f^{(2i-1)}(a)\Big) \\   &&\qquad\qquad\qquad \,+\,c_{m+1}h^{2m+2}(b-a)f^{(2m+2)}(\xi_{m,k}),  \endaligned

gdzie \displaystyle h=(b-a)/k, \displaystyle \xi_{m,k}\in[a,b], a \displaystyle c_i są pewnymi stałymi liczbowymi. Mamy \displaystyle c_1=-1/12, \displaystyle c_2=-1/720 i, ogólnie, \displaystyle c_i=B_i/(2i)!, gdzie \displaystyle B_i są tzw. liczbami Bernoulliego.

Dowód tego lematu pominiemy.

Formułę Eulera-Maclaurina można przepisać w postaci

\displaystyle S(f)\,-\,\bar T_k(f)\,=\,\sum_{i=1}^{m} c^{(0)}_i(f)\,k^{-2i}     \,+\,c^{(0)}_{m+1,k}(f)\,k^{-(2m+2)},

gdzie \displaystyle c^{(0)}_i(f)=c_i(b-a)^{2i}(f^{(2i-1)}(b)-f^{(2i-1)}(a)), \displaystyle 1\le i\le m, oraz \displaystyle c^{(0)}_{m+1,k}(f)=c_{m+1}(b-a)^{2m+2}f^{(2m+2)}(\xi_{m+1,k}). Zauważmy przy tym, że jeśli \displaystyle f\in F^{2m+1}_M([a,b]), to współczynniki \displaystyle c^{(0)}_{m+1,k}(f) są wspólnie ograniczone przez \displaystyle c_{m+1}(b-a)^{2m+2}M.

Definiując teraz kwadraturę

\displaystyle \bar T^1_k(f)\,=\,\frac{4\,\bar T_{2k}(f)\,-\,\bar T_k(f)}{3},

dla \displaystyle f\in C^{(4)}([a,b]) mamy

\displaystyle \aligned S(f)\,-\,\bar T^1_k(f) &= \frac{4\,(S(f)-\bar T_{2k}(f)-    (S(f)-\bar T_k(f))}{3} \\   &= \frac 43\left(\frac{c^{(0)}_1(f)}{4k^2}+       \frac{c^{(0)}_{2,2k}(f)}{4^2k^4}\right)\,-\,       \frac 13\left(\frac{c^{(0)}_1(f)}{k^2}+       \frac{c^{(0)}_{2,k}(f)}{k^4}\right) \\   &= \frac{c^{(1)}_{2,k}(f)}{k^4},  \endaligned

gdzie \displaystyle c^{(1)}_{2,k}(f)=(1/12)c^{(0)}_{2,2k}(f)-(1/3)c^{(0)}_{2,k}(f) i jest wspólnie ograniczone dla \displaystyle f\in F^3_M([a,b]). Kwadratura \displaystyle T^1_k ma więc optymalny w \displaystyle F^3_M([a,b]) rząd zbieżności \displaystyle k^{-4}. Proces ten można kontynuować dalej tworząc kolejne kwadratury o coraz to wyższym rzędzie zbieżności. Dokładniej, połóżmy \displaystyle \bar T^0_k(f)=\bar T_k(f) oraz, dla \displaystyle s\ge 1,

\displaystyle    \bar T^s_k(f)\,=\,\frac   {4^s\,\bar T^{s-1}_{2k}(f)\,-\,\bar T^{s-1}_k(f)}{4^s-1}.

Wtedy, dla \displaystyle f\in F^{2m+1}_M([a,b]), rząd zbieżności kwadratury \displaystyle \bar T^m_k wynosi \displaystyle k^{-(2m+2)}. Rzeczywiście, sprawdziliśmy, że jest to prawdą dla \displaystyle m=0,1. Niech \displaystyle m\ge 2. Postępując indukcyjnie dla \displaystyle s=1,2,\ldots,m mamy

\displaystyle \aligned \lefteqn{  S(f)\,-\,\bar T^s_k(f) \;=\;     \frac{ 4^s(S(f)-\bar T^{s-1}_{2k}(f))-    (S(f)-\bar T^{s-1}_k(f)) }{ 4^s\,-\,1 } } \\   &&=\; \left( 4^s\,\left(\sum_{i=s}^{m}c_i^{(s-1)}(f)(2k)^{-2i}+       c_{m+1,2k}^{(s-1)}(f)(2k)^{-(2m+2)}\right)\right. \\   && \left.\quad \,-\,\left( \sum_{i=s}^mc_i^{(s-1)}(f)k^{-2i}+        c_{m+1,k}^{(s-1)}(f)k^{-(2m+2)}       \right) \right)\,\frac{1}{4^s\,-\,1}\\   &&=\; \sum_{i=s+1}^m c_i^{(s)}(f)k^{-2i}\,+\,         c_{m+1,k}^{(s)}(f)k^{-(2m+2)}, \endaligned

ponieważ współczynniki przy \displaystyle k^{-2s} redukują się. \displaystyle c_i^{(s)}(f) są tutaj pewnymi nowymi stałymi, a \displaystyle c_{m+1,k}^{(s)}(f) może być w klasie \displaystyle F^{2m+1}_M([a,b]) ograniczona przez stałą niezależną od \displaystyle f. Ostatecznie, dla \displaystyle s=m mamy więc

\displaystyle S(f)\,-\,\bar T^m_k(f)\,=\,c_{m+1,k}^{(m)}(f)k^{-(2m+2)}

i w klasie \displaystyle F^{2m+1}_M([a,b])

\displaystyle |S(f)\,-\,\bar T^m_k(f)|\,\le\,c_m\, k^{-(2m+2)}

dla pewnej stałej \displaystyle c_m niezależnej od \displaystyle f.

Zauważmy jeszcze, że \displaystyle \bar T^m_k wykorzystuje \displaystyle n=k2^m+1 wartości \displaystyle f w punktach równoodległych na \displaystyle [a,b], co oznacza, że w terminach \displaystyle n rząd zbieżności wynosi też \displaystyle n^{-(2m+2)}, a więc jest optymalny w klasie \displaystyle F^{2m+1}_M([a,b]).

Kwadratury \displaystyle \bar T^s_k nazywane są kwadraturami Romberga. Dla danej funkcji \displaystyle f można je łatwo konstruować, budując następującą tablicę trójkątną:

\displaystyle    \begin{array} {cccccc}    \bar T^0_1(f) \\    \bar T^0_2(f) &\bar T^1_1(f) \\    \bar T^0_4(f) &\bar T^1_2(f) &\bar T^2_1(f) \\    \bar T^0_8(f) &\bar T^1_4(f) &\bar T^2_2(f)                                  &\bar T^3_1(f) \\    \vdots &\vdots &\vdots &\vdots &\ddots \\    \bar T^0_{2^s}(f) &\bar T^1_{2^{s-1}}(f)           &\bar T^2_{2^{s-2}}(f) &\bar T^3_{2^{s-3}}(f)           &\cdots &\bar T^s_1(f),    \end{array}

której kolumny tworzone są zgodnie z powyższymi wzorami.

Kwadratury adaptacyjne

jak wcześniej zauważyliśmy, błąd kwadratury prostej zależy m.in. od wielkości pochodnej \displaystyle f^{(r+1)} funkcji podcałkowej. Odpowiednia kwadratura złożona wydaje się tego nie zauważać i zagęszcza podział przedziału całkowania jednostajnie, podczas gdy naturalnym i prostym wydaje się pomysł gęstszego podziału tam gdzie \displaystyle |f^{(r+1)}(x)| jest "duża" i rzadszego tam, gdzie \displaystyle f^{(r+1)}(x)| jest "mała". Nasz entuzjazm do tego pomysłu może jednak skutecznie ostudzić uwaga, że algorytm na wejściu zwykle nie dostaje żadnej informacji o \displaystyle f^{(r+1)}. Okazuje się, że mimo wszystko nie stoimy na straconej pozycji. Algorytm obliczający całkę dysponuje na każdym pewną dodatkową informacją o \displaystyle f w postaci jej wartości w pewnych punktach; następny punkt (podział przedziału całkowania) może więc być wybrany na podstawie tych wartości.

Metody uzależniające swoje działanie od konkretnego zadania, które rozwiązują (w naszym przypadku od funkcji podcałkowej) nazywamy ogólnie metodami adaptacyjnymi.

Zauważmy, że poznana wcześniej metoda bisekcji przybliżonego znajdowania zera funkcji jest typową metodą adaptacyjną. Zobaczymy teraz, na przykładzie adaptacyjnej kwadratury Simpsona, jak można wykorzystać adaptację w problemie numerycznego całkowania.

Niech, tak jak poprzednio, \displaystyle \bar P_k będzie złożoną kwadraturą Simpsona z równym podziałem przedziału całkowania na \displaystyle k podprzedziałów, zastosowaną na odcinku \displaystyle [a,b]. W szczególności, \displaystyle \bar P_1=P jest prostą kwadraturą Simpsona. Wtedy

\displaystyle    S(f)-\bar P_2(f)\,=\,-\frac{(b-a)^5}{2280\cdot 2^4}f^{(4)}(\xi_2),

oraz

\displaystyle \aligned \bar P_1(f) - \bar P_2(f) &= (S(f)-\bar P_2(f))\,-\,(S(f)-\bar P_1(f)) \\        &= \frac{(b-a)^5}{2280}\left(f^{(4)}(\xi_1)-\frac 1{16}f^{(4)}(\xi_2)\right), \endaligned

\displaystyle \xi_1,\xi_2\in [a,b]. Załóżmy teraz, że \displaystyle f^{(4)} ma stały znak na \displaystyle [a,b] oraz przedział ten jest na tyle mały, że \displaystyle f^{(4)} jest "prawie stała". Wtedy \displaystyle f^{(4)}(\xi_1)-f^{(4)}(\xi_2)/16\approx 15\cdot f^{(4)}(\xi_2)/16, a stąd otrzymujemy estymator błędu

\displaystyle  S(f) - \bar P_2(f) \approx  -\frac{1}{15}\cdot (\bar P_1(f) -\bar P_2(f)).

Ta przybliżona równość jest podstawą adaptacyjnej kwadratury Simpsona, może bowiem posłużyć do oszacowania błędu na podprzedziałach.

Załóżmy teraz, że chcemy obliczyć wartość całki z dokładnością \displaystyle \varepsilon>0. Obliczamy \displaystyle \bar P_1(f), \displaystyle \bar P_2(f) i sprawdzamy, czy \displaystyle |\bar P_1(f)-\bar P_2(f)|/15\le\varepsilon. Jeśli tak, to \displaystyle \bar P_2(f) jest ostateczną aproksymacją całki na \displaystyle [a,b], a jeśli nie, to dzielimy przedział na dwa podprzedziały \displaystyle [a,(a+b)/2] i \displaystyle [(a+b)/2,b] i powtarzamy procedurę dla obu podprzedziałów z tolerancją błędu \displaystyle \varepsilon/2. Cały proces można zgrabnie zapisać za pomocą funkcji rekurencyjnej.

Algorytm Adaptacyjna kwadratura Simpsona


adaptiveSimpson(a,b,f,e) 
{ 
	P1 = Simpson(a,b,f); 
	P2 = Simpson(a,(a+b)/2,f) + Simpson((a+b)/2,b,f);
	if ( abs(P1-P2) < 15*e) 
		return( P2 );
	else 
		return( adaptiveSimpson(a,(a+b)/2,f,e/2) + adaptiveSimpson((a+b)/2,b,f,e/2) );
}

Zauważmy, że funkcja ta zakończy działanie. Rzeczywiście, na podprzedziale długości \displaystyle h chcemy obliczać całkę z dokładnością \displaystyle \varepsilon h/(b-a), a ponieważ różnica \displaystyle |P1-P2| jest rzędu \displaystyle h^5, kryterium kończenia procedury będzie spełnione dla każdego \displaystyle h dostatecznie małego. Podziały nie mogą więc następować w nieskończoność.

Trochę gorzej sprawa przedstawia się z błędem. Algorytm bazuje bowiem na jego estymatorze. Jeśli po zakończeniu algorytmu mamy podział na podprzedziały \displaystyle a=x_0<x_1<\cdots <x_n=b oraz estymator działa poprawnie na każdym podprzedziale, to błąd można w przybliżeniu oszacować przez

\displaystyle \sum_{j=1}^n\varepsilon\cdot\frac{x_j-x_{j-1}}{b-a}\,=\,\varepsilon.

Z drugiej strony, możemy czasem trafić wyjątkowo "złośliwą" funkcję. Np. jeśli \displaystyle f(a+j(b-a)/4)=0 dla \displaystyle 0\le j\le 4, to już na początku estymator (fałszywie!) twierdzi, że błąd jest zerowy i kwadratura zwróci zero mimo, że rzeczywista wartość całki może znacznie różnić się od zera. Istnieją pewne techniki, które przynajmniej częściowo zapobiegają tego typu zjawiskom, ale nie będziemy ich tutaj omawiać.

Uwarunkowanie całkowania

Zadanie całkowania zadanej funkcji jest, podobnie jak zadanie sumy dwóch liczb (w końcu całkowanie ma wiele wspólnego z sumowaniem!), bardzo dobrze bezwględnie uwarunkowane, natomiast uwarunkowanie względne może być nawet patologicznie duże. Dokładniej,

Twierdzenie O uwarunkowaniu zadania całkowania

Niech \displaystyle f będzie funkcją całkowalną. Wtedy

\displaystyle  \mbox{cond} _{abs}(S,f) = 1

oraz

\displaystyle  \mbox{cond} _{rel}(S,f) = \frac{S(|f|)}{|S(f)|},

gdzie błąd argumentu liczymy w normie \displaystyle ||f|| = \int_a^b|f(x)|\, dx.

Dowód

Biorąc zaburzoną funkcję \displaystyle \tilde{f} taką, że \displaystyle ||\tilde{f} - f||\leq \epsilon mamy

\displaystyle |S(\tilde{f}) - S(f)| = |\int_a^b (\tilde{f} - f)(x) \, dx| \leq \epsilon,

skąd wynika teza.

image:End_of_proof.gif

W szczególności zadanie całkowania będzie źle uwarunkowane względnie, gdy wartość całki jest bliska zeru, ale sama funkcja przyjmuje duże co do modułu wartości.

Biblioteki

W Octave dostępne są jedynie procedury całkujące funkcje skalarne jednej zmiennej na odcinku:

\displaystyle  I = \int_a^b f(x)\, dx.

Robi to funkcja DQAGP ze znakomitego pakietu QUADPACK. Najlepiej od razu posłużmy się przykładem.

Przykład: Prosta całka funkcji jednej zmiennej

Przypuśćmy, że chcemy obliczyć całkę \displaystyle I = \int_0^1 F(x)\, dx, gdzie np. \displaystyle F(x) = \sin(23x) + (1-x^2)^{-1/2}. W tym celu najpierw implementujemy \displaystyle F w Octave:

function y = F(x)
	y = sin(23*x)+1/sqrt(1-x^2);
endfunction

Aby teraz obliczyć całkę \displaystyle I = \int_0^1 F(x)\, dx, wystarczy wywołać

I = quad("F", 0, 1);

W rzeczywistości, podobnie jak w przypadku funkcji fsolve, funkcja quad zwraca więcej informacji, można jej także przekazać dodatkowe parametry. I tak, jeśli chcemy ustawić poziom tolerancji błędu obliczenia całki:

\displaystyle  |I - \text{\lstoct{quad(...)}}| \leq \max\{\text{ATOL}, \text{RTOL}\cdot |I|\}

z wartościami \displaystyle \text{ATOL = 1e-3} i \displaystyle \text{RTOL = 1e-6}, to wywołamy funkcję przekazując jej te parametry następująco:

quad("F", 0, 1, [1e-3, 1e-6]);

Musimy jednak pamiętać, by pojęcia tolerancji "błędu" nie traktować zbyt dosłownie: tym, co naprawdę kontroluje quad podczas wyznaczania wartości całki, jest jedynie pewien estymator błędu, dlatego wartość tolerancji należy zawsze wybierać w sposób konserwatywny, czyli z pewnym zapasem bezpieczeństwa, np.

Jeśli chcemy wyznaczyć wartość całki z błędem bezwzględnym na poziomie \displaystyle 10^{-6}, ustawimy -- na wszelki wypadek -- ATOL = 1e-7, a nie, prostodusznie, ATOL = 1e-6... Musimy także pamiętać, że choć są bardzo mało prawdopodobne do spotkania w praktyce, to jednak istnieją wyuzdane funkcje, dla których estymator błędu może dać całkowicie fałszywe wartości, przez co i obliczona całka może być obarczona dowolnie wielkim błędem.


QUADPACK

Właściwie jedynym klasycznym pakietem, jaki mamy do dyspozycji jest ponaddwudziestoletni QUADPACK. Jest to zestaw kilkunastu procedur fortranowskich, służących obliczaniu typowych całek jednowymiarowych:

Typ całki Procedura QUADPACKa
\displaystyle \int_a^b f(x) DQNG, DQAG, DQAGS, DQAGP
\displaystyle \int_{-\infty}^{\infty} f(x) DQAGI
\displaystyle \int_a^b f(x)\cos(\omega x) DQAWO
\displaystyle \int_a^b \dfrac{f(x)}{x-c} DQAWC

oraz spora liczba podstawowych kwadratur, na których oparto te ogólne. Nazwy procedur rozszyfrowuje się podobnie jak nazwy procedur LAPACKa, zatem

  • przedrostek D w nazwie każdej procedury wymienionej w tabeli (np. DQAGI) oznacza, że będzie działać na liczbach typu double (całkując funkcję \displaystyle f zwracającą wartości tego samego typu). Gdybyśmy chcieli użyć pojedynczej precyzji, użylibyśmy nazwy procedury bez przedrostka.
  • Kolejna litera, Q, oczywiście oznacza kwadraturę (Quadrature).
  • Trzecia litera --- A lub N --- oznacza, odpowiednio, kwadraturę adaptacyjną lub nieadaptacyjną. Jak wiadomo, w praktyce lepiej sprawdzają się kwadratury adaptacyjne, potrafiące w jakiejś mierze dostosować się do przebiegu funkcji podcałkowej. Kwadratury nieadaptacyjne nie mają tej własności. Mogą natomiast, dla pewnych funkcji podcałkowych, okazać się tańsze, warto więc je stosować, gdy wiemy a priori, że adaptacja niewiele pomoże: np. do wolnozmiennych funkcji.
  • Pozostałe litery precyzują typ liczonej całki i zakres ingerencji użytkownika; G --- "zwykła" całka, bez wagi, W --- całka z wagą, O --- dla funkcji silnie oscylujących, C --- wartość główna całki (tzw. całka Cauchy'ego), I --- przedział nieskończony, S --- możliwe osobliwości, P --- użytkownik poda listę punktów, gdzie są osobliwości.

GSL

Biblioteka GSL reimplementuje podstawowe procedury QUADPACKa w języku C. Procedury GSL mają nazwy analogiczne, jak procedury QUADPACKa, ale z przedrostkiem gsl_integration, jak w poniższym przykładzie, gdzie wywołamy odpowiednik procedury DQAG: funkcję gsl_integration_qag.

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_integration.h>
 
double F(double X, void * param) /* wrapper dla funkcji sin(x)/x */
{
	return(sin(X)/X);
}
 
int main(void)
{
	gsl_function f; /* argument z funkcją podcałkową */
 
	double A,ABSERR,B, EPSABS,EPSREL,RESULT;
	int IER,NEVAL;
 
	gsl_integration_workspace *workspace;
 
	int KEY, LIMIT; 
 
	/* przygotowujemy argument z funkcją podcałkową */
	f.function = &F;
	
	A = 0.0E0; B = 10*M_PI; /* przedział całkowania */
	EPSABS = 0.0E0; EPSREL = 1.0E-3; /* tolerancja błędu */
 
	/* parametry specyficzne dla QAG */
	KEY = 1; /* tzn. użyj minimalnej liczby punktów kwadratury bazowej */
	LIMIT = 100; /* maksymalny podział przedziału całkowania */
	workspace  = gsl_integration_workspace_alloc(LIMIT);
	
	/* całkujemy: QAG! */
 
	IER = gsl_integration_qag(&f, A, B, EPSABS, EPSREL, 
					LIMIT, KEY, workspace, &RESULT, &ABSERR);
 
	if (IER != 0)
		fprintf(stderr,"GSL_QAG: Kłopoty z całkowaniem\n");
	fprintf(stderr,"Całka: %g Est. błąd: %g IER: %d\n", RESULT, ABSERR,  IER);
 
	gsl_integration_workspace_free(workspace);
	
	return(0);
}

W powyższym przykładzie specjalnie pozostawiliśmy oznaczenia wykorzystywane w poprzednim programie. Jak widać, funkcje całkujące GSL mają bardzo podobną składnię do odpowiadających im funkcji QUADPACKa.

Miłym rozszerzeniem funkcjonalności jest możliwość przekazywania parametrów do wnętrza funkcji podcałkowej.

Różniczkowanie

<<< Powrót do strony głównej przedmiotu Metody numeryczne

Z zadaniem numerycznego różniczkowania zadanej funkcji spotykamy się często w numeryce. Rzeczywiście, jeśli przypomnimy sobie metodę siecznych, była to po prostu metoda Newtona, w której pochodną przybliżono pewnym ilorazem różnicowym:

\displaystyle x_{k+1} = x_k - \frac{f(x_k)}{g_k},
gdzie
\displaystyle g_k = \frac{f(x_k)-f(x_{k-1})}{x_k - x_{k-1}} \approx f'(x_k).

Zauważmy, że nie jest to jedyny możliwy wzór na przybliżoną metodę Newtona, równie dobrze(? --- to się dopiero okaże!) moglibyśmy wziąć

\displaystyle  g_k = \frac{f(x_k+h)-f(x_k)}{h}

dla dostatecznie małego \displaystyle h.

Podobne formuły są także konieczne do konstrukcji metod numerycznego rozwiązywania równań różniczkowych, gdzie w naturalny sposób pojawia się konieczność operowania pochodną nieznanej funkcji.

Metody różnicowe

Rozważmy najprostszy sposób aproksymacji pochodnej \displaystyle f'(x), oparty na różnicy dzielonej w przód, gdyż ze wzoru Taylora

\displaystyle f(x+h) = f(x) + f'(x)h  + O(h^2),

pomijając człony rzędu \displaystyle h^2, dostajemy przybliżenie

\displaystyle f'(x) \approx \frac{f(x+h)-f(x)}{h},

a dokładniej,

\displaystyle f'(x) =  \frac{f(x+h)-f(x)}{h} + O(h).

Podobną jakość aproksymacji dostaniemy, biorąc różnicę dzieloną w tył,

\displaystyle f'(x) = \frac{f(x)-f(x-h)}{h} + O(h).

Nietrudno przekonać się, że wzięcie średniej arytmetycznej tych dwóch aproksymacji daje tzw. różnicę centralną, która ma wyższy rząd aproksymacji, gdyż

\displaystyle f'(x) =  \frac{f(x+h)-f(x-h)}{2h} + O(h^2),

co znaczy, że dwukrotnie zmniejszając \displaystyle h, powinniśmy się spodziewać aż czterokrotnego zmniejszenia błędu aproksymacji pochodnej!

Jeśli chcemy uzyskać jeszcze wyższy rząd aproksymacji pochodnej, często jako wyrażenie aproksymujące przyjmuje się pochodną wielomianu interpolacyjnego. Rzeczywiście, niech \displaystyle w_n będzie wielomianem interpolującym funkcję \displaystyle f w parami różnych węzłach \displaystyle x_0< \cdots < x_n, tzn.

\displaystyle  w_n(x) = \sum_{i=0}^{n}f(x_i) l_i(x),

gdzie \displaystyle l_i(x)wielomianami bazowymi Lagrange'a. Wtedy

\displaystyle  f'(x) \approx w_n'(x) = \sum_{i=0}^{n}f(x_i) l_i'(x),

przy czym można wykazać, że zachodzi

Twierdzenie O błędzie aproksymacji pochodnej za pomocą pochodnej wielomianu interpolacyjnego

Niech \displaystyle w_n będzie wielomianem interpolującym funkcję \displaystyle f\in C^{n+2}[a,b] w równoodległych węzłach \displaystyle a, a+h, a+2h, \ldots, b, gdzie \displaystyle h = (b-a)/n. Wtedy zachodzi

\displaystyle   f'(x) - w_n'(x) = O(h^n) \quad \forall x\in [a,b]

Wszystkie wprowadzone powyżej metody interpolacji oparte na wielomianie Taylora w rzeczywistości dadzą się sprowadzić do pochodnej wielomianu interpolacyjnego. Rzeczywiście,

  • różnica w przód to aproksymacja \displaystyle f'(x) pochodną wielomianu opartego na węzłach \displaystyle x i \displaystyle x+h,
  • różnica w tył to aproksymacja \displaystyle f'(x) pochodną wielomianu opartego na węzłach \displaystyle x i \displaystyle x-h,
  • różnica centralna to aproksymacja \displaystyle f'(x) pochodną wielomianu opartego na węzłach \displaystyle x-h i \displaystyle x+h; w tym ostatnim przypadku widzimy także, że powyższe twierdzenie nie zawsze jest ostre, bo dla różnicy centralnej byliśmy w stanie uzyskać wyższy niż minimalny gwarantowany przez twierdzenie rząd aproksymacji.

Łatwo także --- korzystając z powyższego --- wyprowadzić nowe wzory na aproksymację pochodnej, przykładowo,

\displaystyle f'(x) = \frac{1}{2h} ( 3f(x) - 4f(x-h) + f(x-2h)) + O(h^2)

korzysta tylko z wartości \displaystyle f na lewo od \displaystyle x (jest to więc różniczkowanie wstecz) i też daje kwadratową aproksymację.

Analogicznie możemy aproksymować pochodne wyższych rzędów, np.

\displaystyle f''(x) = \frac{f(x-h) - 2f(x) + f(x+h)}{h^2} + O(h^2).

Tę formułę możemy znów uzyskać na wiele sposobów:

  • wprost ze wzoru Taylora, raz dla \displaystyle f(x+h), a raz dla \displaystyle f(x-h),
  • jako drugą pochodną wielomianu interpolacyjnego,
  • jako złożenie różnicy dzielonej w przód z różnicą dzieloną w tył (co ma naśladować matematyczną zależność, że \displaystyle f''(x) = (f'(x))'.

Namawiamy czytelnika do sprawdzenia, że faktycznie powyższy wzór można tak wyprowadzić.


Uwarunkowanie różniczkowania

W przeciwieństwie do całkowania, zadanie różniczkowania jest źle postawione ze względu na zaburzenie funkcji, gdy jako dopuszczalne zaburzenia przyjmiemy dowolne funkcje różniczkowalne bliskie danej funkcji w sensie normy jednostajnej. Rzeczywiście, jeśli \displaystyle \tilde{f} jest różniczkowalna, to mimo, że \displaystyle ||\tilde{f} -f||_{C(a,b)} \leq \epsilon, wyrażenie \displaystyle ||\tilde{f}' -f'||_{C(a,b)} nie jest ograniczone. Ten fakt jest źródłem praktycznych kłopotów z numerycznym przybliżaniem pochodnej, gdy np. próbujemy numerycznie różniczkować dane empiryczne: ich błąd często jest funkcją szybkozmienną.

Kłopoty numeryczne z różniczkowaniem

Rozważmy przykładowo różnicę w przód dla \displaystyle h>0. Gdyby arytmetyka której używamy miała nieskończoną precyzję, to oczywiście zachodzi

\displaystyle \frac{f(x+h)-f(x)}{h} - f'(x) = O(h),

i przybliżenie byłoby tym lepsze, im mniejsze byłoby \displaystyle h. Jednak w praktyce tak nie będzie, ze względu na fakt, że działamy w arytmetyce skończonej precyzji:

  • dla małych \displaystyle h, mamy \displaystyle f(x+h) -f(x) \approx 0, a więc zachodzi duże ryzyko utraty cyfr przy odejmowaniu
  • dla małych \displaystyle h, może zdarzyć się, że numerycznie \displaystyle fl_\nu(x+h) = fl_\nu(x) i w konsekwencji \displaystyle fl_\nu(\frac{f(x+h)-f(x)}{h}) = 0.

Można więc postawić sobie pytanie, jak dobrać \displaystyle h na tyle małe, by mieć możliwie dobrą aproksymację \displaystyle f'(x), a jeszcze nie odczuć zgubnych skutków wpływu arytmetyki zmiennoprzecinkowej. Formalnie, możemy pytanie postawić w sposób następujący:

Przypuśćmy, że zamiast \displaystyle f(x) wyznaczane jest \displaystyle \tilde{f}(x) = f(x)+\epsilon_x, przy czym \displaystyle |\epsilon_x| \leq \epsilon. Jak dobrać do \displaystyle \epsilon parametr \displaystyle h w taki sposób, by aproksymacja
\displaystyle \frac{\tilde{f}(x+h)-\tilde{f}(x)}{h} \approx f'(x)

była jak najlepsza? Mamy:

\displaystyle \aligned \left| \frac{\tilde{f}(x+h)-\tilde{f}(x)}{h} - f'(x) \right| &= \left| \frac{1}{h}( (f(x+h)+\epsilon_{x+h}) - (f(x) + \epsilon_x)) - f'(x) \right| \\ &\leq  \left| \frac{ f(x+h) - f(x)}{h} - f'(x)\right| + \frac{2\epsilon}{h} \endaligned

Ponieważ pierwszy człon wyrażenia daje się oszacować (dla dostatecznie regularnej funkcji \displaystyle f) przez \displaystyle C\cdot h, to ostatecznie dostajemy

\displaystyle \left| \frac{\tilde{f}(x+h)-\tilde{f}(x)}{h} - f'(x) \right| \leq C_1(h + \frac{1}{h}).

Wyrażenie po prawej stronie jest minimalizowane dla \displaystyle h = \sqrt{\epsilon} i stąd inżynierska reguła:

Jeśli chcesz używać różnicy w przód, powinieneś wziąć \displaystyle h równe co najmniej \displaystyle \sqrt{\epsilon_{ \mbox{mach} }}.

Literatura

W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj rozdział 7.1 -- 7.6 w

  • D. Kincaid, W. Cheney Analiza numeryczna, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.