Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 63 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 65 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 67 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 69 Notice: Undefined variable: base in /home/misc/mst/public_html/lecture.php on line 36 Matematyka obliczeniowa II – 13. Kwadratury interpolacyjne w wielu wymiarach – MIM UW

Zagadnienia

13. Kwadratury interpolacyjne w wielu wymiarach

13.1. Sformułowanie zadania

Ostatnie trzy wykłady poświęcimy numerycznemu całkowaniu funkcji wielu zmiennych. Dokładniej, dla danej funkcji f:[0,1]^{d}\to\mathbb{R} chcemy obliczyć (przybliżyć) wartość

I_{d}(f)\,=\,\int _{{[0,1]^{d}}}f(\vec{x})\, d\vec{x}\;=\;\underbrace{\int _{0}^{1}\int _{0}^{1}\cdots\int _{0}^{1}}_{d}\, f(x_{1},x_{2},\ldots,x_{d})\, dx_{1}dx_{2}\cdots dx_{d}.

Zakładamy, że powyższa całka istnieje. W ogólniejszym sformułowaniu, chcielibyśmy obliczyć całkę z wagą \omega funkcji f:{\mathbb{R}^{d}}\to\mathbb{R}, która jest postaci

I_{{d,\omega}}(f)\,=\,\int _{{\mathbb{R}^{d}}}\, f(\vec{x})\,\omega(\vec{x})\, d\vec{x}.

Waga \omega jest tutaj nieujemna i całkowalna.

Zauważmy, że ograniczenie się w ostatnim przypadku do \mathbb{R}^{d} nie zmniejsza ogólności, gdyż całkę po dowolnym mierzalnym obszarze D\subseteq\mathbb{R}^{d} można wymodelować przyjmując, że waga \omega(\vec{x})=0 dla wszystkich \vec{x}\notin D.

Zadanie całkowania funkcji wielu zmiennych ma ogromne znaczenie praktyczne i dlatego warto znać skuteczne metody numeryczne jego rozwiązywania.

Przykład 13.1

Wycena obecnej wartości wielu instrumentów finansowych, w tym tzw. opcji, opiera się na założeniu, że przyszłe ceny podlegają losowym zmianom kolejnych odcinkach czasowych. Obecna wartość opcji obliczana jest jako wartość oczekiwana funkcji wypłaty. Odpowiada to obliczaniu całki oznaczonej funkcji d zmiennych, gdzie d jest liczbą odciników czasowych. Jest to często całka ze standardową d wymiarową wagą gaussowską postaci

(2\pi)^{{-d/2}}\int _{{\mathbb{R}^{d}}}f(\xi _{1},\ldots,\xi _{d})\exp\Big(-\frac{1}{2}(\xi _{1}^{2}+\cdots,\xi _{d}^{2})\Big)\, d\xi _{1}\cdots d\xi _{d},

przy czym f jest (zwykle skomplikowaną) funkcją wypłaty na końcu okresu, a \xi _{j} reprezentują czynniki losowe w kolejnych odcinkach czasu. Wymiar d może wynosić nawet kilka tysięcy.

Z podstawowego wykładu z metod numerycznych każdy z nas wie jak numerycznie całkować funkcje jednej zmiennej. Stosowane metody w znakomitej większości przypadków sprowadzają się do scałkowania funkcji, która jest kawałkami wielomianem określonego stopnia interpolującym funkcję podcałkową. Pomysł ten może być uogólniony na przypadek funkcji wielu zmiennych. Aby jednak mówić o kwadraturach interpolacyjnych w wielu wymiarach, musimy najpierw zastanowić się nad rozwiązywalnością odpowiedniego zadania interpolacyjnego.

13.2. Interpolacja na siatkach regularnych

13.2.1. Postać wielomianu interpolacyjnego

Niech

a\le t_{1}<t_{2}<\cdots<t_{r}\le b.

Jeśli f jest funkcją jednej zmiennej, f:[a,b]\to\mathbb{R}, to wielomian

p_{f}(x)=\sum _{{j=1}}^{r}f(t_{j})\, l_{j}(x),

gdzie l_{j} jest j-tym wielomianem Lagrange'a,

l_{j}(x)=\prod _{{j\ne i=1}}^{r}\frac{x-t_{j}}{t_{i}-t_{j}},\qquad 1\le j\le r,

(przy czym l_{1}\equiv 1 dla r=1) jest stopnia co najwyżej (r-1) i interpoluje f w punktach t_{j}, tzn. przyjmuje w tych punktach te same wartości co f. W przypadku d\ge 2 możemy podobnie zdefiniować ,,wielowymiarowe” wielomiany Lagrange'a.

W tym celu zakładamy, że na każdej współrzędnej dany jest przedział, a w nim układ r punktów

a^{{(k)}}\le t_{1}^{{(k)}}<t_{2}^{{(k)}}<\ldots<t_{r}^{{(k)}}\le b^{{(k)}},\qquad 1\le k\le d.

Oznaczając przez l_{j}^{{(k)}} odpowiednie wielomiany Lagrange'a jednej zmiennej dla k-tego podziału, definiujemy wielomiany Lagrange'a d zmiennych jako

l_{{j_{1},\ldots,j_{d}}}(x_{1},\ldots,x_{d})=l^{{(1)}}_{{j_{1}}}(x_{1})\, l^{{(2)}}_{{j_{2}}}(x_{2})\cdots l^{{(d)}}_{{j_{d}}}(x_{d})

dla wszystkich 1\le j_{k}\le r, 1\le k\le d. Dla skrócenia zapisu, będziemy dalej używać zapisu wektorowego \vec{j}=(j_{1},\ldots,j_{d}), a 1\le\vec{j}\le d będzie oznaczać, że nierówności zachodzą dla każdej współrzędnej j_{k}, 1\le k\le d. Podobnie, t_{{\vec{j}}}=(t^{{(1)}}_{{j_{1}}},t^{{(2)}}_{{j_{2}}},\ldots,t^{{(d)}}_{{j_{d}}}).

Wielomiany l_{{\vec{j}}} należą do przestrzeni {\mathcal{P}}_{d}^{r} wielomianów d zmiennych postaci

p(\vec{x})=p(x_{1},\ldots,x_{d})=\sum _{{0\le\vec{i}\le r-1}}a_{{\vec{i}}}\cdot x_{1}^{{i_{1}}}x_{2}^{{i_{2}}}\cdots x_{d}^{{i_{d}}},

gdzie a_{{\vec{i}}} są dowolnymi wsoółczynnikami rzeczywistymi. Zauważmy, że p\in{\mathcal{P}}_{d}^{r} wtedy i tylko wtedy gdy p jest wielomianem stopnia co najwyżej (r-1) ze względu na każdą zmienną x_{k}.

Lemat 13.1

Jeśli wielomian p\in{\mathcal{P}}_{d}^{r} zeruje się we wszystkich r^{d} punktach t_{{\vec{j}}}, 1\le\vec{j}\le r, to p jest wielomianem zerowym.

Dowód

Dowód przeprowadzimy przez indukcję ze względu na wymiar d. Dla d=1 lemat jest oczywiście prawdziwy, bo na podstawie zasadniczego twierdzenia algebry niezerowy wielomian stopnia co najwyżej (r-1) nie może mieć r różnych zer.

Niech d\ge 2. Niech a_{{\vec{j}}} będą współczynnikami wielomianu p. Dla ustalonej k zdefiniujmy wielomian p_{k}\in{\mathcal{P}}_{{d-1}}^{r} jako

p_{k}(x_{1},\ldots,x_{{d-1}})\,=\, p(x_{1},\ldots,x_{{d-1}},t^{{(d)}}_{k}).

Wielomian ten zeruje się w r(d-1) punktach t_{{i_{1},\ldots,i_{{d-1}}}}. Zapisując go w postaci

p_{k}(x_{1},\ldots,x_{{d-1}})\,=\,\sum _{{0\le i_{1},\ldots,i_{{d-1}}\le d-1}}b^{{(k)}}_{{i_{1},\ldots,i_{{d-1}}}}\cdot x_{1}^{{i_{1}}}\cdots x_{{d-1}}^{{i_{{d-1}}}},

gdzie współczynniki

b^{{(k)}}_{{i_{1},\ldots,i_{{d-1}}}}\,=\,\sum _{{0\le i_{d}\le d-1}}a_{{\vec{i}}}\cdot(t_{k}^{{(d)}})^{{i_{d}}},

oraz stosując założenie indukcyjne mamy, że b^{{(k)}}_{{i_{1},\ldots,i_{{d-1}}}}=0. A więc dla wszystkich wyborów indeksów i_{1},\ldots,i_{{d-1}} wielomian jednej zmiennej \sum _{{i_{d}=0}}^{{d-1}}a_{{\vec{i}}}\cdot t^{{i_{d}}} zeruje się w r punktach t=t^{{(d)}}_{s}. To zaś wymusza a_{{\vec{i}}}=0 dla wszystkich 0\le\vec{i}\le d-1 i w konsekwencji p\equiv 0.

Lemat 13.1 wykorzystamy do pokazania następującego twierdzenia.

Twierdzenie 13.1

Wielomiany l_{{\vec{j}}}, 1\le\vec{j}\le r, tworzą bazę przestrzeni {\mathcal{P}}_{d}^{r}. W szczególności, {\rm dim}({\mathcal{P}}_{d}^{r})=r^{d}.

Dowód

Zauważmy, że podobnie jak w przypadku d=1,

l_{{\vec{j}}}(t_{{\vec{i}}})\,=\,\left\{\begin{array}[]{ll}1,&\mbox{jeśli}\;\vec{i}=\vec{j},\\
0,&\mbox{w przeciwnym przypadku}.\end{array}\right.

Stąd, jeśli kombinacja liniowa p=\sum _{{\vec{j}}}\alpha _{{\vec{j}}}l_{{\vec{j}}} jest wielomianem zerowym to dla wszystkich \vec{i}

0=p(t_{{\vec{i}}})=\sum _{{\vec{j}}}\alpha _{{\vec{j}}}l_{{\vec{j}}}\,(t_{{\vec{i}}})=\alpha _{{\vec{i}}},

czyli układ \{ l_{{\vec{i}}}:\, 1\le\vec{i}\le r\} jest liniowo niezależny. Z drugiej strony, układ ten rozpina {\mathcal{P}}_{d}^{r}, bo dla dowolnego wielomianu p z tej przestrzeni mamy

p=\sum _{{1\le\vec{j}\le r}}p(t_{{\vec{j}}})\, l_{{\vec{j}}}. (13.1)

Rzeczywiście, w przeciwnym przypadku różnica wielomianu p i prawej strony (13.1) byłaby niezerowym wielomianem w {\mathcal{P}}_{d}^{r}, który zeruje się we wszystkich r^{d} punktach t_{{\vec{i}}}. To zaś przczyłoby lematowi 13.1.

Stąd już jeden krok do następującego wniosku podsumowującego nasze dotychczasowe rozważania. Niech

D\,=\,[a^{{(1)}},b^{{(1)}}]\times[a^{{(2)}},b^{{(2)}}]\times\cdots\times[a^{{(d)}},b^{{(d)}}]

będzie d wymiarowym prostokątem.

Wniosek 13.1

Dla dowolnej funkcji f:D\to\mathbb{R} wielomian

p_{f}(\vec{x})\,=\,\sum _{{0\le\vec{j}\le r}}f(t_{{\vec{j}}})\, l_{{\vec{j}}}(\vec{j})

jest jedynym wielomianem w {\mathcal{P}}_{d}^{r} interpolującym f w punktach t_{{\vec{j}}} tzn. takim, że

p_{f}(t_{{\vec{j}}})=f(t_{{\vec{j}}})

dla wszystkich 1\le\vec{j}\le r.

13.2.2. Błąd interpolacji

Zastanówmy się teraz jaki jest błąd otrzymanej interpolacji. Dla uproszczenia będziemy od teraz zakładać, że D jest kostką d wymiarową, tzn. wszystkie krawędzie mają tą samą długość, którą oznaczymy przez H, a węzły na każdej współrzędnej

t_{j}^{{(k)}}=a^{{(k)}}+u_{j}H,\qquad 1\le j\le r,

gdzie

0\le u_{1}<u_{2}<\cdots<u_{r}\le 1

jest pewną ustaloną siatką na odcinku jednostkowym.

W przypadku skalarnym, o ile funkcja f jest r-krotnie różniczkowalna w sposób ciągły, to

f(x)-p_{f}(x)\,=\,(x-t_{1})(x-t_{2})\cdots(x-t_{r})\frac{f^{{(r)}}(\xi)}{r!},

przy czym \xi\in[a,b] zależy od x. Stąd w szczególności mamy

|f(x)-p_{f}(x)|\,\le\,\frac{(b-a)^{r}}{r!}\| f^{{(r)}}\| _{\infty}, (13.2)

gdzie \| f^{{(r)}}\| _{\infty}=\max _{{a\le t\le b}}|f^{{(r)}}(t)|. Aby wyprowadzić formułę na bład interpolacji w przypadku wielowymiarowym, będziemy potrzebować pewnego prostego uogólnienia ostatniego wzoru.

Załóżmy, że zamiast dokładnych wartości f(t_{i}) mamy jedynie wartości przybliżone y_{i} takie, że błąd

|y_{i}-f(t_{i})|\,\le\,\delta,\qquad 1\le i\le r. (13.3)

Niech dalej \tilde{p}_{f} będzie wielomianem stopnia co najwyżej (r-1) interpolującym dane przybliżone y_{i} w punktach t_{i}. Ponieważ (p_{f}-\tilde{p}_{f}) jest wielomianem interpolującym dane f(t_{j})-y_{j}, na podstawie wzoru (13.1) mamy

|p_{f}(x)-\tilde{p}_{f}(x)|\,\le\,\delta\cdot\sum _{{i=1}}^{r}|l_{i}(x)|\,\le\,\delta\cdot S_{r},

gdzie S_{1}=1, a dla r\ge 2

S_{r}=\max _{{0\le z\le 1}}\sum _{{i=1}}^{r}\prod _{{i\ne j=1}}^{r}\left|\frac{z-u_{j}}{u_{i}-u_{j}}\right|.

Stąd i z formuły na błąd interpolacji dla dokładnych danych otrzymujemy

|f(x)-\tilde{p}_{f}(x)|\,\le\,|f(x)-p_{f}(x)|+|p_{f}(x)-\tilde{p}_{f}(x)|\,\le\,\frac{(b-a)^{r}}{r!}\| f^{{(r)}}\| _{\infty}\,+\,\delta\cdot S_{r}. (13.4)

Wprowadzimy jeszcze klasę {\mathcal{F}}_{r}(D) funkcji f:D\to\mathbb{R}, które w całej swojej dziedzinie są r-krotnie różniczkowalne w sposób ciągły ze względu na każdą zmienną. Dla f\in{\mathcal{F}}_{r}(D) definiujemy

B_{{r}}(f)\,=\,\max _{{1\le i\le d}}\left\{\,\left\|\frac{\partial^{r}f}{\partial x_{1}^{r}}\right\| _{\infty},\ldots,\left\|\frac{\partial^{r}f}{\partial x_{d}^{r}}\right\| _{\infty}\,\right\}.
Twierdzenie 13.2

Niech D=[a^{{(1)}},a^{{(1)}}+H]\times\cdots\times[a^{{(d)}},a^{{(d)}}+H]. Jeśli f\in{\mathcal{F}}_{r}(D) to dla każdego \vec{x}\in D błąd interpolacji

|f(\vec{x})-p_{f}(\vec{x})|\,\le\,\frac{H^{r}}{r!}\, C_{{r,d}}\, B_{r}(f),

gdzie C_{{1,d}}=d, a dla r\ge 2

C_{{r,d}}\,=\,\frac{S_{r}^{d}-1}{S_{r}-1}.
Dowód

Rozpatrzymy tylko r\ge 2 pozostawiając przypadek r=1 jako proste ćwiczenie.

Dla d=1 nierówność w tezie jest równoważna (13.2). Załóżmy więc, że d\ge 2. Ponieważ dla każdego ustalonego t^{{(d)}}_{k} wielomian (d-1) zmiennych p_{f}(x_{1},\ldots,x_{{d-1}},t^{{(d)}}_{k}) jest wielomianem interpolacyjnym dla funkcji (d-1) zmiennych f(x_{1},\ldots,x_{{d-1}},t^{{(d)}}_{k}), na podstawie założenia indukcyjnego mamy

|f(x_{1},\ldots,x_{{d-1}},t^{{(d)}}_{k})-p_{f}(x_{1},\ldots,x_{{d-1}},t^{{(d)}}_{k})|\,\le\,\frac{H^{r}}{r!}\, B_{{r}}(f)\left(\frac{S_{r}^{{d-1}}-1}{S_{r}-1}\right). (13.5)

Zauważmy, że dla ustalonych z kolei pierwszych (d-1) współrzędnych x_{1},\ldots,x_{{d-1}} wielomian p_{f}(x_{1},\ldots,x_{{d-1}},t) jest wielomianem jednej zmiennej t interpoluącym funkcję jednej zmiennej f(x_{1},\ldots,x_{{d-1}},t) w punktach t^{{(d)}}_{k} na podstawie danych zaburzonych na poziomie \delta równym prawej stronie (13.5). Stąd i z (13.4) ostatecznie otrzymujemy

\displaystyle|f(\vec{x})-p_{f}(\vec{x})| \displaystyle\le \displaystyle\frac{H^{r}}{r!}\, B_{r}(f)\,+\,\delta\cdot S_{r}
\displaystyle= \displaystyle\frac{H^{r}}{r!}\, B_{r}(f)\,\left(1+\frac{S_{r}^{{d-1}}-1}{S_{r}-1}\, S_{r}\right)
\displaystyle= \displaystyle\frac{H^{r}}{r!}\, B_{r}(f)\,\left(\frac{S_{r}^{d}-1}{S_{r}-1}\right).

13.3. Kwadratury interpolacyjne

13.3.1. Kwadratury proste

Jesteśmy już dobrze uzbrojeni w mechanizm interpolacyjny i możemy zdefiniować wielowymiarowe kwadratury interpolacyjne dla całkowania funkcji f:D\to\mathbb{R} zdefiniowanych na kostce

D=[a^{{(1)}},a^{{(1)}}+H]\times\cdots\times[a^{{(d)}},a^{{(d)}}+H].

Kwadratury te dane są równością

Q_{{r,d}}(f)\,=\,\int _{{D}}p_{f}(\vec{x})\, d\vec{x}, (13.6)

gdzie p_{f}\in{\mathcal{P}}_{d}^{r} jest wielomianem interpolującym f w punktach t_{{\vec{j}}}, 1\le\vec{j}\le r.

Chociaż postać (13.6) kwadratury znakomicie nadaje się do rozważań teoretycznych, nie jest jednak praktyczna ze względu na obliczenia. Zauważmy, że

\displaystyle Q_{{r,d}}(f) \displaystyle= \displaystyle\int _{{D}}\sum _{{\vec{j}}}f(t_{{\vec{j}}})l_{{\vec{j}}}(\vec{x})\, d\vec{x}\;=\;\sum _{{\vec{j}}}f(t_{{\vec{j}}})\,\int _{{D}}l_{{\vec{j}}}(\vec{x})\, d\vec{x}
\displaystyle= \displaystyle H^{d}\cdot\sum _{{\vec{j}}}f(t_{{\vec{j}}})\,\prod _{{k=1}}^{d}\left(\int _{0}^{1}l_{{j_{k}}}(u)\, du\right),

gdzie l_{j} jest j-tym wielomianem Lagrange'a dla punktów u_{1},u_{2},\ldots,u_{d}. Stąd, wprowadzając oznaczenie

a_{k}\,=\,\int _{0}^{1}l_{k}(u)\, du,

kwadraturę interpolacyjną można zapisać w postaci

Q_{{r,d}}(f)\,=\, H^{d}\cdot\sum _{{1\le j_{1},\ldots,j_{d}\le r}}a_{{j_{1}}}a_{{j_{2}}}\cdots a_{{j_{d}}}\cdot f(t^{{(1)}}_{{j_{1}}},t^{{(2)}}_{{j_{2}}},\ldots,t^{{(d)}}_{{j_{d}}}).

Zauważmy, że a_{k} są współczynnikami jednowymiarowej kwadratury interpolacyjnej Q_{r}(f)=\sum _{{k=1}}^{r}a_{k}f(t_{k}) opartej na punktach u_{k}, przybliżającej całkę \int _{0}^{1}f(u)\, du. Mówiąc inaczej, zdefiniowana przez nas wielowymiarowa kwadratura interpolacyjna jest d-produktem tensorowym wybranej kwadratury jednowymiarowej.

Na koniec tego podrozdziału podamy oszacowanie błędu kwadratury Q_{{r,d}}. Ponieważ

\int _{{D}}f(\vec{x})\, d\vec{x}\,-\, Q_{{r,d}}(f)\;=\;\int _{{D}}\big(f(\vec{x})-p_{f}(\vec{x})\big)\, d\vec{x},

z twierdzenia 13.2 natychmiast otrzymujemy następujący wniosek.

Wniosek 13.2

Jeśli f\in{\mathcal{F}}_{r}(D) to błąd kwadratury interpolacyjnej Q_{{r,d}} jest ograniczony przez

\left|\int _{{D}}f(\vec{x})\, d\vec{x}\,-\, Q_{{r,d}}(f)\right|\;\le\;\frac{H^{{r+d}}}{r!}\, C_{{r,d}}\, B_{r}(f).

13.3.2. Kwadratury złożone

Podobnie jak w przypadku funkcji jednej zmiennej, definiujemy kwadratury złożone dla funkcji wielu zmiennych. Dla uproszczenia zakładamy, że całkujemy po kostce jednostkowej [0,1]^{d}.

Dla danego n wprowadzamy podział kostki na n^{d} podkostek

\left[\frac{i_{1}-1}{n},\frac{i_{1}}{n}\right]\times\left[\frac{i_{2}-1}{n},\frac{i_{2}}{n}\right]\times\cdots\left[\frac{i_{d}-1}{n},\frac{i_{d}}{n}\right],\qquad 1\le i_{k}\le n,\quad 1\le k\le d.

Następnie na każdej podkostce stosujemy prostą kwadraturę interpolacyjną opartą na siatce regularnej składającej się z r^{d} punktów. Skonstruowaną w ten sposób kwadraturę złożoną oznaczymy przez Q^{{(n)}}_{{r,d}}.

Przykład 13.2

Jeśli bazową kwadraturą jednowymiarową jest reguła punktu środkowego,

Q_{1}(f)\,=\,(b-a)\cdot f\left(\frac{a+b}{2}\right),

to wynikową kwadraturą złożoną na [0,1]^{d} jest po prostu reguła prostokątów

Q^{{(n)}}_{{r,d}}(f)\,=\,\left(\frac{1}{n}\right)^{d}\cdot\sum _{{1\le i_{1},,\ldots,i_{d}\le n}}f\left(\frac{i_{1}-1/2}{n},\cdots,\frac{i_{d}-1/2}{n}\right).

Nasze rozważania wieńczy twierdzenie o błędzie kwadratury złożonej, które natychmiast wynika z wniosku 13.2 oraz sposobu konstrukcji kwadratury.

Twierdzenie 13.3

Kwadratura złożona Q^{{(n)}}_{{r,d}} korzysta z co najwyżej

N\,=\,(r\, n)^{d}

wartości funkcji f. Jeśli f\in{\mathcal{F}}_{r}([0,1]^{d}) to jej błąd

\left|\int _{{[0,1]^{d}}}f(\vec{x})\, d\vec{x}\,-\, Q^{{(n)}}_{{r,d}}(f)\right|\,\le\,\left(\frac{1}{N}\right)^{{r/d}}\left(\frac{r^{r}}{r!}\right)\, C_{{r,d}}\, B_{r}(f).

13.4. Przekleństwo wymiaru

Złożone kwadratury interpolacyjne mogą być z powodzeniem stosowane dla niskich wymiarów, powiedzmy d=2,3. Dla dużych wymiarów d mają one bowiem tą niepożądaną własność, że liczba węzłów rośnie wykładniczo szybko wraz z zagęszczaniem siatki. Nawet jeśli weźmiemy po 2 punkty na każdej współrzędnej to całkowita liczba punktów siatki regularnej wyniesie 2^{d}. Pamiętamy, że w wielu praktycznych zastosowaniach d może sięgać nawet kilku tysięcy. W takich przypadkach obliczenie wartości kwadratury jest zadaniem praktycznie niewykonalnym.

To jednak nie koniec złych wiadomości. Przyjrzyjmy się jeszcze błędowi złożonej kwadratury interpolacyjnej. Twierdzenie 13.3 mówi, że błąd ten jest ograniczony z góry proporcjonalnie do N^{{-r/d}}, gdzie N jest liczbą wszystkich użytych punktów. To drugi powód do niepokoju, uzasadniony poniższym przykładem.

Przykład 13.3

Załóżmy, że chcemy całkować funkcję 360 zmiennych i jako kwadraturę bazową stosujemy kwadraturę Simpsona, dla której r=4. Górne ograniczenie błędu sugeruje, że aby być pewnym wyniku z dokładnością 10^{{-2}} to musimy obliczyć wartości funkcji w aż 10^{{180}} punktach. Czy naprawdę jest aż tak źle?

Rzeczywiście jest tak źle, a nawet gorzej. Okazuje się, że rzędu zbieżości N^{{-r/d}} nie da się poprawić w klasie funkcji {\mathcal{F}}_{r}([0,1]^{d}). Mówi o tym następujące twierdzenie.

Twierdzenie 13.4

Istnieje c=c_{{r,d}}>0 o następującej własności. Dla dowolnej aproksymacji całki wykorzystującej N wartości funkcji można znaleźć funkcję f\in{\mathcal{F}}_{r}([0,1]^{d}) dla której B_{r}(f)=1, a błąd aproksymacji całki wynosi co najmniej c\, N^{{-r/d}}.

Dowód

Załóżmy, że dana aproksymacja całki oblicza wartości funkcji w punktach \vec{t}_{j}, 1\le j\le N. Dowód twierdzenia polega na konstrukcji dwóch funkcji, f_{+} i f_{-}, które zerują się we wszystkich \vec{t}_{j} (a tym samym ich całki są aproksymowane tą samą liczbą), dla których B_{r}(f_{+})=1=B_{r}(f_{-}), ale różnica całek

\int _{{[0,1]^{d}}}(f_{+}-f_{-})(\vec{t})\, d\vec{t}\,\ge\, 2c\, N^{{-r/d}},

dla pewnej c niezależnej od f i d. Wtedy, przynajmniej dla jednej z tych funkcji błąd aproksymacji całki wynosi co najmniej cN^{{-r/d}}.

Wybierzmy n taką, że n^{d}>N i skonstruujmy na [0,1]^{d} regularną siatkę składającą się z n^{d} kostek, każda o krawędzi długości h=1/n.

Niech dalej \phi:\mathbb{R}\to\mathbb{R} będzie dowolną funkcją r-krotnie różniczkowalną w sposób ciągły spełniającą następujące warunki:

  1. \phi(x)=0 dla x\notin(0,1),

  2. \phi^{{(j)}}(0)=0=\phi^{{(j)}}(1) dla 0\le j\le r,

  3. \int _{0}^{1}\phi(t)\, dt=:a>0.

Każdej kostce

K_{{\vec{i}}}:=[(i_{1}-1)h,i_{1}h]\times\cdots\times[(i_{d}-1)h,i_{d}h]

naszej regularnej siatki przyporządkujemy funkcję

\phi _{{\vec{i}}}(x_{1},\ldots,x_{d}):=h^{r}\prod _{{k=1}}^{d}\phi(x_{k}/h-i_{k}).

Zauważmy, że B_{r}(\phi _{{\vec{i}}})=1 oraz

\int _{{K_{{\vec{i}}}}}\phi _{{\vec{i}}}(\vec{t})\, d\vec{t}\,=\, a^{d}\, h^{{r+d}}.

Jasne jest, że istnieje co najmniej n^{d}-N multi-indeksów \vec{i} (kostek) takich, że żaden z punktów \vec{t}_{j} nie należy do wnętrza K_{{\vec{i}}}. Oznaczmy zbiór takich indeksów przez S i zdefiniujmy funkcje

f_{+}:=\sum _{{\vec{i}\in S}}\phi _{{\vec{i}}},\qquad f_{-}:=-f_{+}.

Wtedy obie funkcje zerują się w \vec{t}_{j}, B_{r}(f_{+})=1=B_{r}(f_{-}), oraz

\int _{{[0,1]^{d}}}f_{+}(\vec{t})\, d\vec{t}\,=\,-\int _{{[0,1]^{d}}}f_{-}(\vec{t})\, d\vec{t}\,\ge\,(n^{d}-N)\, a^{d}\, h^{{r+d}}.

Podstawiając n=\lceil N^{{1/d}}(1+d/r)^{{1/d}}\rceil dostajemy ostatecznie

\int _{{[0,1]^{d}}}f_{+}(\vec{t})\, d\vec{t}\,\ge\, c\, N^{{-r/d}},\qquad\mbox{gdzie}\quad c=\frac{da^{d}}{r2^{{r+d}}(1+d/r)^{{1+r/d}}}.

Opisane zjawisko nosi nazwę przekleństwa wymiaru.

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

Projekt współfinansowany przez Unię Europejską w ramach Europejskiego Funduszu Społecznego.

Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.