MN14

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

$S(f)=\underset{a}{\overset{b}{\int }}f(x)dx,$

gdzie $-\mathrm{\infty }<a<b<+\mathrm{\infty }$, a $f$ należy do pewnej klasy $F$ funkcji rzeczywistych określonych i całkowalnych w sensie Riemanna na całym przedziale $[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 $f$, a w niektórych przypadkach również jej pochodnych, o ile istnieją. Dokładna całka $S(f)$ będzie więc w ogólności przybliżana wartością $A(f)$, która zależy tylko od wartości $f$ i ew. jej pochodnych w skończonej liczbie punktów.

Kwadratury

Kwadraturami nazywamy funkcjonały liniowe Parser nie mógł rozpoznać (nieznana funkcja „\toR”): {\displaystyle \displaystyle Q:F\toR} postaci

$Q(f)=\sum _{i=0}^n{a}_{i}f(x_i),$

albo ogólniej

$Q(f)=\sum _{i=0}^k\sum _{j=0}^{n_i-1}{a}_{i,j}{f}^{(j)}(x_i),$

gdzie $x_i$ są punktami z $[a,b]$, a $a_i$ (albo $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 $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 $f$, otrzymana w ten sposób wartość też będzie zależeć tylko od tej informacji, a w konsekwencji funkcjonał wynikowy będzie postaci ja wyżej. Są to tzw. kwadratury interpolacyjne.

Definicja

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

$Q^I(f)=\underset{a}{\overset{b}{\int }}{w}_f(x)dx,$

gdzie $w_f$ jest wielomianem interpolacyjnym funkcji $f$ stopnia co najwyżej $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 $l_i$ (zob. (Uzupelnic: Lagrbaza )), otrzymujemy

$Q^I(f)=\underset{a}{\overset{b}{\int }}\sum _{i=0}^{n}f(x_i)l_i(x)dx=\sum _{i=0}^{n}f(x_i)\underset{a}{\overset{b}{\int }}{l}_i(x)dx,$

a stąd i z postaci $l_i$,

$a_i=\underset{a}{\overset{b}{\int }}\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,$

$0\le i\le n$.

Podamy teraz kilka przykładów.

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

$Q_0^I(f)=(b-a)f\left(\frac{a+b}{2}\right).$

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

$Q_1^I(f)=T(f)=\frac{b-a}{2}(f(a)+f(b)).$

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

$Q_2^I(f)=P(f)=\frac{b-a}{6}(f(a)+4f(\frac{a+b}{2})+f(b)).$

Zauważmy, że kwadratury trapezów i parabol są oparte na węzłach jednokrotnych i równoodległych, przy czym $x_0=a$ i $x_n=b$. Ogólnie, kwadratury interpolacyjne oparte na węzłach równoodległych Parser nie mógł rozpoznać (błąd składni): {\displaystyle \displaystyle x_i=a+(b-a)i/ń} , $0\le i\le n$, nazywamy kwadraturami Newtona--Cotesa.

Błąd kwadratur interpolacyjnych

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

Twierdzenie

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

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

W klasie $F_M^n([a,b])$, maksymalny błąd kwadratury $Q^I$ wynosi

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

Dowód\quad Korzystając ze znanego nam już wzoru na błąd interpolacji wielomianowej z Lematu Uzupelnic: leblint , mamy

$S(f)-Q^I(f)=\underset{a}{\overset{b}{\int }}(x-x_0)(x-x_1)\cdots (x-x_n)f(x_0,x_1,\dots ,x_n,x)dx.$

Stąd, jeśli $f\in F_M^n([a,b])$ to

$|S(f)-Q^I(f)|\le \underset{a}{\overset{b}{\int }}(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 $F_M^n([a,b])$ wynika bezpośrednio z oszacowania (Uzupelnic: blkwad ). Aby pokazać ograniczenie dolne zauważmy, że dla funkcji $g$ takiej, że $g^{(n+1)}$ przyjmuje na przedziałach $(a,x_0)$, $(x_0,x_1)$, $\dots$, $(x_n,b)$ naprzemiennie wartości $M$ i $-M$ mamy

$|S(g)-Q^I(g)|=\frac{M}{(n+1)!}\underset{a}{\overset{b}{\int }}|(x-x_0)(x-x_1)\cdots (x-x_n)|dx.$

Co prawda, $g$ nie jest w $F_M^n[a,b]$, ale może być dla dowolnego $ϵ>0$ przybliżana funkcjami $f_{ϵ}\in F_M^n([a,b])$ w ten sposób, że całka

$\underset{a}{\overset{b}{\int }}|(x-x_0)\cdots (x-x_n)(f-g{)}^{(n+1)}(x)|dx\le ϵ(n+1)!.$

Zapisując $f_{ϵ}=g+(f_{ϵ}-g)$ mamy

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \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 $ϵ$ daje dowód twierdzenia. $◻$

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

Twierdzenie

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

$S(f)-T(f)=-\frac{(b-a{)}^3}{12}{f}^{(2)}({\xi }_1).$

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

$S(f)-P(f)=-\frac{(b-a{)}^5}{2280}{f}^{(4)}({\xi }_2).$

(${\xi }_1,{\xi }_2\in [a,b]$).

Dowód\quad (i) Ze wzoru (Uzupelnic: blkwad )

$S(f)-T(f)=\underset{a}{\overset{b}{\int }}(x-a)(x-b)f(a,b,x)dx.$

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

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \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 $c,{\xi }_1\in [a,b]$.

(ii) Niech $w_{f,2}\in {\mathrm{\Pi }}_2$ i $w_{f,3}\in {\mathrm{\Pi }}_3$ będą wielomianami interpolacyjnymi funkcji $f$ odpowiednio dla węzłów $a,b,(a+b)//2$ oraz $a,b,(a+b)//2,(a+b)//2$. Wtedy

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

Wobec

$\underset{a}{\overset{b}{\int }}(x-a)\left(\frac{a+b}{2}\right)(x-b)dx=0$

mamy

$P(f)=\underset{a}{\overset{b}{\int }}{w}_{f,2}(x)dx=\underset{a}{\overset{b}{\int }}{w}_{f,3}(x)dx.$

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

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \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 $(x-a)(x-(a+b)//2{)}^2(x-b)$ jest niedodatni na $[a,b]$, możemy znów zastosować twierdzenie o wartości średniej. Mamy

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \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. $◻$

Kwadratury złożone

Podobnie jak w przypadku zadania interpolacji 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 $f$.

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

$\overline{Q}(f)=\sum _{i=0}^n(t_{i+1}-t_i)f(x_i),$

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

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

$\overline{Q}(f)=\underset{a}{\overset{b}{\int }}{\overline{w}}_f(x)dx,$

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

Twierdzenie

Błąd kwadratury złożonej $\overline{Q}(f)$ w klasie $F_M^r([a,b])$ jest ograniczony przez

$\underset{f\in F_M^r([a,b])}{\sup }|S(f)-\overline{Q}(f)|\le \frac{(b-a{)}^{r+2}}{k^{r+1}}\frac{M}{(r+1)!}\le C(\frac{1}{n}{)}^{r+1},$

gdzie

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

Dowód\quad Twierdzenie to jest bezpośrednim wnioskiem z Twierdzenia Uzupelnic: twblkw . Mamy bowiem

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \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. $◻$

W klasie $F_M^r([a,b])$, błąd kwadratur złożonych jest rzędu $n^{-(r+1)}$, czyli jest taki sam jak błąd interpolacji kawałkami wielomianowej. I tak jak przy interpolacji można pokazać, że błąd każdej innej metody całkowania korzystającej jedynie z wartości funkcji w $n$ punktach nie może w klasie $F_M^r([a,b])$ maleć szybciej niż $n^{-(r+1)}$, zob. U. Uzupelnic: optzlo . 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 ${\overline{T}}_k$ i parabol ${\overline{P}}_k$. Powstają one przez zastosowanie na każdym przedziale $[t_{i-1},t_i]$ odpowiednio kwadratur trapezów $T$ i parabol $P$. Jak łatwo się przekonać,

${\overline{T}}_k(f)=\frac{b-a}{k}(\frac{f(a)+f(b)}{2}+\sum _{j=1}^{k-1}f(\frac{j}{k}\left)\right),$

oraz

${\overline{P}}_k(f)=\frac{b-a}{3k}(\frac{f(a)+f(b)}{2}+\sum _{j=1}^{k-1}f(\frac{j}{k})+2\sum _{j=1}^{k}f(\frac{2j-1}{2k}\left)\right).$

Twierdzenie

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

$S(f)-{\overline{T}}_k(f)=-\frac{(b-a{)}^2}{12k^2}{f}^{(2)}({\xi }_1).$

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

$S(f)-{\overline{P}}_k(f)=-\frac{(b-a{)}^4}{2280k^4}{f}^{(4)}({\xi }_2).$

Dowód

Dla kwadratury trapezów mamy

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \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)^2}{12 k^2}\sum_{i=1}^k \frac{b-a}k f^{(2)}(\alpha_i) \,=\, -\frac{(b-a)^2}{12 k^2} f^{(2)}(\xi_1), \endaligned}

a dla kwadratury parabol podobnie

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \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)^4}{2280 k^4}\sum_{i=1}^k \frac{b-a}k f^{(4)}(\beta_i) \,=\, -\frac{(b-a)^4}{2280 k^4} f^{(4)}(\xi_2). \endaligned}

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

Przyspieszanie zbieżności kwadratur

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

${\overline{T}}_{2k}=\frac{1}{2}({\overline{T}}_k(f)+\frac{b-a}{k}\sum _{i=1}^{k}f(\frac{2i-1}{2k}\left)\right).$

Pozwala on obliczyć ${\overline{T}}_{2k}(f)$ na podstawie ${\overline{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ń ${\overline{T}}_{2^s}(f)$ ($s\ge 0$) całki $S(f)$. Jest to szczególnie istotne wtedy, gdy nie mamy żadnej informacji a priori o $\Vert f^{″}{\Vert }_{C([a,b])}$, a przez to nie potrafimy oszacować liczby $n$ węzłów, dla której osiągniemy pożądaną dokładność, zob. U. Uzupelnic: kwas .

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ż $n^{-2}$. Okazuje się jednak, że kwadratury ${\overline{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ę

Colin Maclaurin
Zobacz biografię

Lemat Formuła Eulera-Maclaurina

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

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \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 $h=(b-a)//k$, ${\xi }_{m,k}\in [a,b]$, a $c_i$ są pewnymi stałymi liczbowymi. Mamy $c_1=-1//12$, $c_2=-1//720$ i, ogólnie, $c_i=B_i//(2i)!$, gdzie $B_i$ są tzw. liczbami Bernoulliego. $◻$

Dowód tego lematu pominiemy.

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

$S(f)-{\overline{T}}_k(f)=\sum _{i=1}^m{c}_i^{(0)}(f)k^{-2i}+c_{m+1,k}^{(0)}(f)k^{-(2m+2)},$

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

Definiując teraz kwadraturę

${\overline{T}}_k^1(f)=\frac{4{\overline{T}}_{2k}(f)-{\overline{T}}_k(f)}{3},$

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

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \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 $c_{2,k}^{(1)}(f)=(1//12)c_{2,2k}^{(0)}(f)-(1//3)c_{2,k}^{(0)}(f)$ i jest wspólnie ograniczone dla $f\in F_M^3([a,b])$. Kwadratura $T_k^1$ ma więc optymalny w $F_M^3([a,b])$ rząd zbieżności $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 ${\overline{T}}_k^0(f)={\overline{T}}_k(f)$ oraz, dla $s\ge 1$,

${\overline{T}}_k^s(f)=\frac{4^s{\overline{T}}_{2k}^{s-1}(f)-{\overline{T}}_k^{s-1}(f)}{4^s-1}.$

Wtedy, dla $f\in F_M^{2m+1}([a,b])$, rząd zbieżności kwadratury ${\overline{T}}_k^m$ wynosi $k^{-(2m+2)}$. Rzeczywiście, sprawdziliśmy, że jest to prawdą dla $m=0,1$. Niech $m\ge 2$. Postępując indukcyjnie ze względu na $s=1,2,\dots ,m$ mamy

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \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 $k^{-2s}$ redukują się. $c_i^{(s)}(f)$ są tutaj pewnymi nowymi stałymi, a $c_{m+1,k}^{(s)}(f)$ może być w klasie $F_M^{2m+1}([a,b])$ ograniczona przez stałą niezależną od $f$. Ostatecznie, dla $s=m$ mamy więc

$S(f)-{\overline{T}}_k^m(f)=c_{m+1,k}^{(m)}(f)k^{-(2m+2)}$

i w klasie $F_M^{2m+1}([a,b])$

$|S(f)-{\overline{T}}_k^m(f)|\le c_m{k}^{-(2m+2)}$

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

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

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

$\begin{array}{c}{\overline{T}}_1^0(f)\\ {\overline{T}}_2^0(f)& {\overline{T}}_1^1(f)\\ {\overline{T}}_4^0(f)& {\overline{T}}_2^1(f)& {\overline{T}}_1^2(f)\\ {\overline{T}}_8^0(f)& {\overline{T}}_4^1(f)& {\overline{T}}_2^2(f)& {\overline{T}}_1^3(f)\\ ⋮& ⋮& ⋮& ⋮& \ddots \\ {\overline{T}}_{2^s}^0(f)& {\overline{T}}_{2^{s-1}}^1(f)& {\overline{T}}_{2^{s-2}}^2(f)& {\overline{T}}_{2^{s-3}}^3(f)& \cdots & {\overline{T}}_1^s(f),\end{array}$

w której pierwsza kolumna jest tworzona indukcyjnie zgodnie ze wzorem (Uzupelnic: indtrap ), a kolejne zgodnie z (Uzupelnic: indRom ).

Biblioteki

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

$I=\underset{a}{\overset{b}{\int }}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ę $I=\underset{0}{\overset{1}{\int }}F(x)dx$, gdzie np. $F(x)=\sin (23x)+(1-x^2{)}^{-1/2}$. W tym celu najpierw implementujemy $F$ w Octave:

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

Aby teraz obliczyć całkę $I=\underset{0}{\overset{1}{\int }}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:

Parser nie mógł rozpoznać (błąd składni): {\displaystyle \displaystyle |I - \text{\lstoct{quad(...)}}| \leq \max\{\text{ATOL}, \text{RTOL}\cdot |I|\} }

z wartościami $\text{ATOL = 1e-3}$ i $\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 $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.

Dodajmy, że quad doskonale radzi sobie także z funkcjami z osobliwościami (o ile tylko ją o nich uprzedzimy). Przykładowo, scałkujmy funkcję ze złośliwą nieciągłością w zerze:

Parser nie mógł rozpoznać (nieznana funkcja „\begincases”): {\displaystyle \displaystyle f(x) = \begincases \frac{1}{x}, \qquad x\neq 0,\\ 10^6, \qquad x=0. \endcases }

Oczywiście, $\underset{-1}{\overset{1}{\int }}f(x)dx=0$, tymczasem definiując

function y = osobliwa(x)
if (x != 0.0)
	y = 1/x;
else
	y = 1e6;
endif
endfunction

quad("osobliwa", -1, 1)

daje w rezultacie

ABNORMAL RETURN FROM DQAGP ans = -81.098

-- błędny wynik i (na szczęście!) także komunikat o jakichś problemach procedury całkującej. Jeśli jednak wspomożemy quad informacją o tym, że w $x=0$ czai się osobliwość, wszystko przebiegnie już gładko:

quad("osobliwa", -1, 1, [1e-8 1e-8], 0) ans = 0

(1e-8 jest żądaną tolerancją błędu, równą wartości domyślnej).

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, w tym kilka ogólnego stosowania:

Typ całki	Procedura QUADPACKa
$\underset{a}{\overset{b}{\int }}f(x)$	DQNG, DQAG, DQAGS, DQAGP
$\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}f(x)$	DQAGI
$\underset{a}{\overset{b}{\int }}f(x)\cos (\omega x)$	DQAWO
$\underset{a}{\overset{b}{\int }}\frac{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ę $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, natomiast są tańsze, warto więc je stosować w szczególnych przypadkach, gdy wiemy a priori, że adaptacja niewiele pomoże: np. do wolno zmiennych 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 --- oscylacyjną, 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, zob. {GSL-reference}, gdzie zostały one bardzo przystępnie i szczegółowo opisane. 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, to była to po prostu metoda Newtona, w której pochodną przybliżono pewnym ilorazem różnicowym:

$x_{k+1}=x_k-\frac{f(x_k)}{g_k},$

gdzie

$g_k=\frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}}\approx f^{\prime }(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ąć

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

dla dostatecznie małego $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 $f^{\prime }(x)$, oparty na różnicy dzielonej w przód, gdyż ze wzoru Taylora

$f(x+h)=f(x)+f^{\prime }(x)h+O(h^2),$

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

$f^{\prime }(x)\approx \frac{f(x+h)-f(x)}{h},$

a dokładniej,

$f^{\prime }(x)=\frac{f(x+h)-f(x)}{h}+O(h).$

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

$f^{\prime }(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ż

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

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

Metody interpolacyjne

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 $w_n$ będzie wielomianem interpolującym funkcję $f$ w parami różnych węzłach $x_0<\cdots <x_n$, tzn.

$w_n(x)=\sum _{i=0}^{n}f(x_i)l_i(x),$

gdzie $l_i(x)$ są \link{sec:interpol}{wielomianami bazowymi Lagrange'a}. Wtedy

$f^{\prime }(x)\approx {w_n}^{\prime }(x)=\sum _{i=0}^{n}f(x_i){l_i}^{\prime }(x),$

przy czym można wykazać, że

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

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

$f^{\prime }(x)-{w_n}^{\prime }(x)=O(h^n)\mathrm{\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 $f^{\prime }(x)$ pochodną wielomianu opartego na węzłach $x$ i $x+h$,
• różnica w tył to aproksymacja $f^{\prime }(x)$ pochodną wielomianu opartego na węzłach $x$ i $x-h$,
• różnica centralna to aproksymacja $f^{\prime }(x)$ pochodną wielomianu opartego na węzłach $x-h$ i $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ę, przykładowo,

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

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

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

$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 $f(x+h)$, a raz dla $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 $f^{″}(x)=(f^{\prime }(x){)}^{\prime }$.

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

Kłopoty numeryczne

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

$\frac{f(x+h)-f(x)}{h}-f^{\prime }(x)=O(h),$

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

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

Można więc postawić sobie pytanie, jak dobrać $h$ na tyle małe, by mieć możliwie dobrą aproksymację $f^{\prime }(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

$f(x)$

wyznaczane jest

$\tilde{f}(x)=f(x)+{ϵ}_x$

, przy czym

$|{ϵ}_x|\le ϵ$

. Jak dobrać do

$ϵ$

parametr

$h$

w taki sposób, by aproksymacja

$\frac{\tilde{f}(x+h)-\tilde{f}(x)}{h}\approx f^{\prime }(x)$

była jak najlepsza?

Mamy:

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \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 $f$) przez $C\cdot h$, to ostatecznie dostajemy

$|\frac{\tilde{f}(x+h)-\tilde{f}(x)}{h}-f^{\prime }(x)|\le C_1(h+\frac{1}{h}).$

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

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