MN14

Całkowanie

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]$.

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.

\noindent 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).$

\noindent 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)).$

\noindent 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.$

\noindent 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

(i) 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).$

(ii) 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]$).

\noindent 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]$.

\noindent (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)!}.$

\noindent 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ę

Plik:MacLaurin.jpg

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:

\begintable [ht]

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

\caption{Tabela najważniejszych procedur QUADPACKa i ich zastosowań.} \endtable

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 Uzupelnic: table:quadpack (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: 
	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.