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 ewentualnie 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 takiej postaci, jak 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, 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, są natomiast 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.

Numeryczne równania różniczkowe

Dla rozwiązywania równań różniczkowych zwyczajnych,

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned y'(t) &= f(y(t),t), \qquad t\in (t_0,T]\\ y(t_0) &= y_0, \endaligned}

korzystamy w Octave z funkcji lsode. Jak już zdążyliśmy się przyzwyczaić, także ta funkcja jest jedynie zgrabnym interfejsem dla procedury DLSODE z pakietu ODEPACK autorstwa A.Hindmarscha. Domyślnie wybierana jest metoda BDF (dobra dla równań trudniejszych, tzw. sztywnych równań różniczkowych, lecz względnie kosztowna). Należy dodać, że solver kontroluje zachowanie się rozwiązania i sam w trakcie obliczeń może zmienić wartości pewnych swoich parametrów w zależności od bieżących potrzeb.

Wprowadzenie do numerycznych metod rozwiązywania RRZ

LSODE jest pakietem, którego jądrem obliczeniowym są tzw. schematy różnicowe. Nie miejsce tu na dokładną analizę numeryczną metod rozwiązywania równań różniczkowych zwyczajnych (zainteresowanych odsyłam do doskonałej monografii Hairer et. al , bądź do książki Palczewskiego), ograniczymy się tutaj do pewnych intuicji. Będą nam one potrzebne, by móc świadomie wpływać na pewne parametry metody w celu uzyskania lepszego (lub szybszego) rozwiązania (zazwyczaj te dwie własności stoją w opozycji względem siebie).

Jawne schematy różnicowe

Najprostszym schematem aproksymacji równania różniczkowego zwyczajnego jest schemat Eulera: gdy pochodną przybliżymy ilorazem różnicowym, dostajemy (przybliżone) równanie

Leonhard Euler
Zobacz biografię

$\frac{y^h(t+h)-y^h(t)}{h}=f(y^h(t),t).$

(dla podkreślenia, że spełniająca je funkcja jest jedynie przybliżeniem dokładnego rozwiązania $y$, oznaczyliśmy je $y^h$).

Oznaczając dalej $t_1=t_0+h$, $t_2=t_1+h$, ... $t_n=t_{n-1}+h$ zauważamy, że
z powyższej równości wynika wzór rekurencyjny

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned y_{n+1} - y_{n} = h\cdot f_n, \endaligned}

gdzie $y_n=y^h(t_n)\approx y(t_n)$ oraz $f_n=f(y_n,t_n)$. Naturalnie, warunkiem początkowym pozostaje $y_0$. Jest to schemat jawny, gdyż $y_{n+1}$ daje się wyznaczyć wprost z równania -- jest zadane w postaci nieuwikłanej.

Uzasadnieniem dla schematu Eulera jest to, że iloraz $\frac{y_{n+1}-y_n}{h}$ dość dobrze przybliża pochodną, gdy $h$ jest małe:

$\frac{y_{n+1}-y_n}{h}=y^{\prime }(t_n)+O(h).$

Oczywiście pochodną $y^{\prime }(t_n)$ możemy przybliżać na różne sposoby, korzystając np. z wielomianu interpolacyjnego dla $y$. Takie podejście prowadzi do klasy metod wielokrokowych, ogólnej postaci

$\sum _{j=0}^k{\alpha }_k{y}_{n+k}=h\sum _{j=0}^k{\beta }_k{f}_{n+k},$

(metoda $k$-krokowa). Inną ważną klasą metod są metody jednokrokowe tzw. Rungego--Kutty.

Niejawne schematy różnicowe

Bardzo podobnie wyglądający niejawny schemat Eulera,

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned y_{n+1} - y_{n} = h\cdot f_{n+1}, \endaligned}

ma znacząco różne własności. Pierwsza, która rzuca się w oczy to konieczność rozwiązywania równania nieliniowego na kazdym kroku schematu (co znacząco podnosi jego koszt), bo tym razem $y_{n+1}$ zadaje się wzorem niejawnym: $y_{n+1}-y_n=h\cdot f(y_{n+1},t_{n+1})$. Za chwilę pokażemy, że czasem opłaca się zapłacić tę cenę, gdyż schematy niejawne mają pewne bardzo pożądane własności.

Lokalna dokładność schematu

Każdy schemat rozwiązywania

równania zwyczajnego jest schematem przybliżonym, a jakość lokalnego przybliżenia równania schematem można ściśle określić przy użyciu wielkości całkowitej zwanej rzędem schematu: im wyższy, tym większa może być wartość $h$, przy której schemat dostatecznie dobrze przybliżać będzie równanie. Przykładowo, schemat Eulera stosunkowo mało dokładnie aproksymuje równanie dla niezbyt małych $h$: aby uzyskać rozwiązanie obarczone możliwie małym błędem, należy zazwyczaj wybierać dostatecznie mały krok schematu $h$, a to znaczy, że w końcowym efekcie będziemy musieli wykonać bardzo wiele kroków schematu, co podniesie jego całkowity koszt.

Należy pamiętać, że schematy wysokiego rzędu nie muszą prowadzić do uzyskania dokładniejszego rozwiązania. Więcej szczegółów na ten temat Czytelnik odnajdzie w fachowym monografiach, pamiętając, że metody użyte w Octave potrafią do pewnego stopnia automatycznie dopasować krok całkowania i rząd schematu tak, by uzyskać możliwie dobre przybliżenie prawdziwego rozwiązania równania.

Równania sztywne

Schematy jawne okazują się niewystarczająco skuteczne dla pewnych ważnych klas równań, tzw. równań sztywnych. Niektórzy autorzy wręcz definiują równania sztywne jako te, dla których schematy jawne nie są skuteczne tzn. np. schemat jawny z adaptacyjnym wyborem długości kroku musi stosować patologicznie krótki krok.

Rozważmy sprawę sztywności na przykładzie układu równań różniczkowych liniowych

$\frac{dY}{dt}=-\mathrm{\Lambda }Y$

(gdzie $Y(t)=(y_1(t),\dots ,y_K(t){)}^T\in R^K$, a $\mathrm{\Lambda }$ jest macierzą $K\times K$).

Przypuśćmy, że macierz $\mathrm{\Lambda }$ jest diagonalizowalna, a liczby ${\lambda }_1,\dots ,{\lambda }_K$ są jej wartościami własnymi. Gdy wszystkie ${\lambda }_i$ mają ujemne części rzeczywiste, składowe $y_i$ rozwiązania gasną do zera, gdy $t\to \mathrm{\infty }$. Rozsądne jest, aby od schematu wymagać podobnej właściwości, tzn. by składowe $y_{i,n}^h$ rozwiązania przybliżonego miały własność $y_{i,n}^h\to 0$, gdy $n\to \mathrm{\infty }$. Okazuje się, że dla jawnej metody Eulera tak będzie, gdy

$h\le \frac{2}{\underset{i}{\max }{\lambda }_i}.$

A więc długość kroku czasowego jawnej metody Eulera nie może być zbyt wielka, co jest szokujące, bo dla dużych wartości $t$ rozwiązania stabilizują się wokół zera i wydawać by się mogło, że to pozwoliłoby zwiększyć wówczas krok schematu. Niestety tak nie jest i krok schematu musimy dostosować do najmniejszej skali czasowej, która objawia się jedynie na początku ewolucji układu. I to jest cecha charakterystyczna dla ogólnie nazywanych równań sztywnych, nie tylko liniowych jak w powyższym przykładzie.

Tymczasem, dla niejawnego schematu Eulera nie ma w omawianym przykładzie żadnego ograniczenia na krok schematu! Ogólniej, schematy niejawne, (czy to wielokrokowe typu Adamsa--Bashfortha lub BDF, czy niejawne metody Rungego--Kutty) lepiej nadają się do rozwiązywania równań sztywnych, niż odpowiadające im analogiczne schematy jawne. Ceną, jaką płacimy za stosowanie schematu niejawnego, a tym samym za uzyskanie mniejszych ograniczeń stabilnościowych na długość użytego kroku schematu, jest większy (być może nawet znacznie większy) koszt wykonania jednego kroku schematu, bo tym razem na każdym kroku schematu musimy rozwiązywać równanie algebraiczne, w ogólnym przypadku nieliniowe. lsode domyślnie używa właśnie schematu BDF, by użytkownik był z góry zabezpieczony na wypadek sztywności rozwiązywanego zadania.

Niestety równania sztywne pojawiają się dość często w zastosowaniach, co generalnie wiąże się z obecnością różnych skal czasowych ewolucji poszczególnych składników rozważanego modelu. Sztywne są m.in. układy równań różniczkowych zwyczajnych powstałe po dyskretyzacji (przestrzennej) równań różniczkowych cząstkowych typu parabolicznego. Również wiele układów równań różniczkowych zwyczajnych spotykanych w zastosowaniach, np. w modelach reakcji chemicznych, jest z natury sztywnych. Podkreślmy raz jeszcze, że z analizy teoretycznej i z praktycznych doświadczeń wynika, że dla tej klasy zadań niezbędne jest stosowanie schematów niejawnych, np. typu BDF, w których $y_{n+1}$ będzie zadane w sposób uwikłany.

Biblioteki

Działanie instrukcji Octave'a lsode najlepiej prześledzić na konkretnym przykładzie.

Rozwiązujemy równanie Van der Pola

Rozważmy równanie różniczkowe zwyczajne rzędu drugiego, opisujące tak zwany oscylator Van der Pola:

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned y'' + \epsilon (y^2-1)y' + y = 0. \endaligned}

Ponieważ metody dla równań zwyczajnych zazwyczaj pisze się dla równań pierwszego rzędu (i takie są procedury Octave), zapiszemy najpierw to równanie w postaci układu równań pierwszego rzędu:

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned y_1' &= y_2\\ y_2' &= - \epsilon (y_1^2 - 1) y_2 - y_1. \endaligned}

Jest to równanie różniczkowe postaci

$y^{\prime }(t)=f(y(t),t),t\in [t_0,T],$

gdzie u nas $y$ jest wektorem o dwóch współrzędnych, $y=(y_1,y_2{)}^T$. Aby znaleźć numeryczne rozwiązanie tego równania w Octave, najpierw zdefiniujemy M-funkcję, która oblicza prawą stronę $f(y,t)$:

 
function ydot = vanderpol(y,t)
global epsilon;

ydot = [y(2);
        epsilon*(1-y(1)^2)*y(2) - y(1)];
endfunction

Zwróćmy uwagę na wykorzystanie zmiennej globalnej epsilon, dzięki której będzie możliwe przekazanie wartości parametru $ϵ$ do wnętrza funkcji vanderpol(y,t). Procedura rozwiązująca równanie wymaga funkcji, której argumentami są tylko i wyłącznie liczba $t$ i wektor $y$ o tylu współrzędnych ile jest równań w układzie.

Choć trudno w to uwierzyć, to jednak w zasadzie najtrudniejszy krok mamy za sobą! Teraz, aby rozwiązać równanie, wystarczy wywołać funkcję

lsode("vanderpol",y0,t),

gdzie y0 będzie wektorem wartości początkowych rozwiązania: $y(t_0)=y_0$, a wektor t - zestawem kolejnych punktów (pierwszym punktem w zestawie musi być oczywiście chwila początkowa $t_0$) w przedziale czasowym $[t_0,T]$, w których zechcemy poznać wartości rozwiązania: naturalnie, aby dobrze przedstawić rozwiązanie graficznie, musimy wyznaczyć je w dostatecznie wielu punktach.

CVODE

Z dobrych bibliotek napisanych w C należy wymienić pakiet CVODE. Omówienie jego możliwości wykracza jednak poza zakres niniejszego wykładu.