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,1dR chcemy obliczyć (przybliżyć) wartość

Idf=0,1dfxdx=010101dfx1,x2,,xddx1dx2dxd.

Zakładamy, że powyższa całka istnieje. W ogólniejszym sformułowaniu, chcielibyśmy obliczyć całkę z wagą ω funkcji f:RdR, która jest postaci

Id,ωf=Rdfxωxdx.

Waga ω jest tutaj nieujemna i całkowalna.

Zauważmy, że ograniczenie się w ostatnim przypadku do Rd nie zmniejsza ogólności, gdyż całkę po dowolnym mierzalnym obszarze DRd można wymodelować przyjmując, że waga ωx=0 dla wszystkich xD.

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π-d/2Rdfξ1,,ξdexp-12ξ12+,ξd2dξ1dξd,

przy czym f jest (zwykle skomplikowaną) funkcją wypłaty na końcu okresu, a ξ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

at1<t2<<trb.

Jeśli f jest funkcją jednej zmiennej, f:a,bR, to wielomian

pfx=j=1rftjljx,

gdzie lj jest j-tym wielomianem Lagrange'a,

ljx=ji=1rx-tjti-tj,1jr,

(przy czym l11 dla r=1) jest stopnia co najwyżej r-1 i interpoluje f w punktach tj, tzn. przyjmuje w tych punktach te same wartości co f. W przypadku d2 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

akt1k<t2k<<trkbk,1kd.

Oznaczając przez ljk odpowiednie wielomiany Lagrange'a jednej zmiennej dla k-tego podziału, definiujemy wielomiany Lagrange'a d zmiennych jako

lj1,,jdx1,,xd=lj11x1lj22x2ljddxd

dla wszystkich 1jkr, 1kd. Dla skrócenia zapisu, będziemy dalej używać zapisu wektorowego j=j1,,jd, a 1jd będzie oznaczać, że nierówności zachodzą dla każdej współrzędnej jk, 1kd. Podobnie, tj=tj11,tj22,,tjdd.

Wielomiany lj należą do przestrzeni Pdr wielomianów d zmiennych postaci

px=px1,,xd=0ir-1aix1i1x2i2xdid,

gdzie ai są dowolnymi wsoółczynnikami rzeczywistymi. Zauważmy, że pPdr wtedy i tylko wtedy gdy p jest wielomianem stopnia co najwyżej r-1 ze względu na każdą zmienną xk.

Lemat 13.1

Jeśli wielomian pPdr zeruje się we wszystkich rd punktach tj, 1jr, 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 d2. Niech aj będą współczynnikami wielomianu p. Dla ustalonej k zdefiniujmy wielomian pkPd-1r jako

pkx1,,xd-1=px1,,xd-1,tkd.

Wielomian ten zeruje się w rd-1 punktach ti1,,id-1. Zapisując go w postaci

pkx1,,xd-1=0i1,,id-1d-1bi1,,id-1kx1i1xd-1id-1,

gdzie współczynniki

bi1,,id-1k=0idd-1aitkdid,

oraz stosując założenie indukcyjne mamy, że bi1,,id-1k=0. A więc dla wszystkich wyborów indeksów i1,,id-1 wielomian jednej zmiennej id=0d-1aitid zeruje się w r punktach t=tsd. To zaś wymusza ai=0 dla wszystkich 0id-1 i w konsekwencji p0.

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

Twierdzenie 13.1

Wielomiany lj, 1jr, tworzą bazę przestrzeni Pdr. W szczególności, dimPdr=rd.

Dowód

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

ljti=1,jeślii=j,0,w przeciwnym przypadku.

Stąd, jeśli kombinacja liniowa p=jαjlj jest wielomianem zerowym to dla wszystkich i

0=pti=jαjljti=αi,

czyli układ li: 1ir jest liniowo niezależny. Z drugiej strony, układ ten rozpina Pdr, bo dla dowolnego wielomianu p z tej przestrzeni mamy

p=1jrptjlj.(13.1)

Rzeczywiście, w przeciwnym przypadku różnica wielomianu p i prawej strony (13.1) byłaby niezerowym wielomianem w Pdr, który zeruje się we wszystkich rd punktach ti. 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=×a1,b1×a2,b2××ad,bd

będzie d wymiarowym prostokątem.

Wniosek 13.1

Dla dowolnej funkcji f:DR wielomian

pfx=0jrftjljj

jest jedynym wielomianem w Pdr interpolującym f w punktach tj tzn. takim, że

pftj=ftj

dla wszystkich 1jr.

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

tjk=ak+ujH,1jr,

gdzie

0u1<u2<<ur1

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

fx-pfx=x-t1x-t2x-trfrξr!,

przy czym ξa,b zależy od x. Stąd w szczególności mamy

fx-pfxb-arr!fr,(13.2)

gdzie fr=maxatbfrt. 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 fti mamy jedynie wartości przybliżone yi takie, że błąd

yi-ftiδ,1ir.(13.3)

Niech dalej p~f będzie wielomianem stopnia co najwyżej r-1 interpolującym dane przybliżone yi w punktach ti. Ponieważ pf-p~f jest wielomianem interpolującym dane ftj-yj, na podstawie wzoru (13.1) mamy

pfx-p~fxδi=1rlixδSr,

gdzie S1=1, a dla r2

Sr=max0z1i=1rij=1rz-ujui-uj.

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

fx-p~fxfx-pfx+pfx-p~fxb-arr!fr+δSr.(13.4)

Wprowadzimy jeszcze klasę FrD funkcji f:DR, 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 fFrD definiujemy

Brf=max1idrfx1r,,rfxdr.
Twierdzenie 13.2

Niech D=a1,a1+H××ad,ad+H. Jeśli fFrD to dla każdego xD błąd interpolacji

fx-pfxHrr!Cr,dBrf,

gdzie C1,d=d, a dla r2

Cr,d=Srd-1Sr-1.
Dowód

Rozpatrzymy tylko r2 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 d2. Ponieważ dla każdego ustalonego tkd wielomian d-1 zmiennych pfx1,,xd-1,tkd jest wielomianem interpolacyjnym dla funkcji d-1 zmiennych fx1,,xd-1,tkd, na podstawie założenia indukcyjnego mamy

fx1,,xd-1,tkd-pfx1,,xd-1,tkdHrr!BrfSrd-1-1Sr-1.(13.5)

Zauważmy, że dla ustalonych z kolei pierwszych d-1 współrzędnych x1,,xd-1 wielomian pfx1,,xd-1,t jest wielomianem jednej zmiennej t interpoluącym funkcję jednej zmiennej fx1,,xd-1,t w punktach tkd na podstawie danych zaburzonych na poziomie δ równym prawej stronie (13.5). Stąd i z (13.4) ostatecznie otrzymujemy

fx-pfxHrr!Brf+δSr
=Hrr!Brf1+Srd-1-1Sr-1Sr
=Hrr!BrfSrd-1Sr-1.

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:DR zdefiniowanych na kostce

D=a1,a1+H××ad,ad+H.

Kwadratury te dane są równością

Qr,df=Dpfxdx,(13.6)

gdzie pfPdr jest wielomianem interpolującym f w punktach tj, 1jr.

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

Qr,df=Djftjljxdx=jftjDljxdx
=Hdjftjk=1d01ljkudu,

gdzie lj jest j-tym wielomianem Lagrange'a dla punktów u1,u2,,ud. Stąd, wprowadzając oznaczenie

ak=01lkudu,

kwadraturę interpolacyjną można zapisać w postaci

Qr,df=Hd1j1,,jdraj1aj2ajdftj11,tj22,,tjdd.

Zauważmy, że ak są współczynnikami jednowymiarowej kwadratury interpolacyjnej Qrf=k=1rakftk opartej na punktach uk, przybliżającej całkę 01fudu. 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 Qr,d. Ponieważ

Dfxdx-Qr,df=Dfx-pfxdx,

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

Wniosek 13.2

Jeśli fFrD to błąd kwadratury interpolacyjnej Qr,d jest ograniczony przez

Dfxdx-Qr,dfHr+dr!Cr,dBrf.

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,1d.

Dla danego n wprowadzamy podział kostki na nd podkostek

i1-1n,i1ni2-1n,i2nid-1n,idn,1ikn,1kd.

Następnie na każdej podkostce stosujemy prostą kwadraturę interpolacyjną opartą na siatce regularnej składającej się z rd punktów. Skonstruowaną w ten sposób kwadraturę złożoną oznaczymy przez Qr,dn.

Przykład 13.2

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

Q1f=b-afa+b2,

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

Qr,dnf=1nd1i1,,,idnfi1-1/2n,,id-1/2n.

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 Qr,dn korzysta z co najwyżej

N=rnd

wartości funkcji f. Jeśli fFr0,1d to jej błąd

0,1dfxdx-Qr,dnf1Nr/drrr!Cr,dBrf.

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 2d. 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ż 10180 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 Fr0,1d. Mówi o tym następujące twierdzenie.

Twierdzenie 13.4

Istnieje c=cr,d>0 o następującej własności. Dla dowolnej aproksymacji całki wykorzystującej N wartości funkcji można znaleźć funkcję fFr0,1d dla której Brf=1, a błąd aproksymacji całki wynosi co najmniej cN-r/d.

Dowód

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

0,1df+-f-tdt 2cN-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 nd>N i skonstruujmy na 0,1d regularną siatkę składającą się z nd kostek, każda o krawędzi długości h=1/n.

Niech dalej ϕ:RR będzie dowolną funkcją r-krotnie różniczkowalną w sposób ciągły spełniającą następujące warunki:

  1. ϕx=0 dla x0,1,

  2. ϕj0=0=ϕj1 dla 0jr,

  3. 01ϕ(t)dt=:a>0.

Każdej kostce

Ki:=i1-1h,i1h××id-1h,idh

naszej regularnej siatki przyporządkujemy funkcję

ϕix1,,xd:=hrk=1dϕxk/h-ik.

Zauważmy, że Brϕi=1 oraz

Kiϕitdt=adhr+d.

Jasne jest, że istnieje co najmniej nd-N multi-indeksów i (kostek) takich, że żaden z punktów tj nie należy do wnętrza Ki. Oznaczmy zbiór takich indeksów przez S i zdefiniujmy funkcje

f+:=iSϕi,f-:=-f+.

Wtedy obie funkcje zerują się w tj, Brf+=1=Brf-, oraz

0,1df+tdt=-0,1df-tdtnd-Nadhr+d.

Podstawiając n=N1/d1+d/r1/d dostajemy ostatecznie

0,1df+tdtcN-r/d,gdziec=dadr2r+d1+d/r1+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.