Kwadratury interpolacyjne

Zbiór wzorów, definicji i najczęściej poruszanych problemów z Analizy.
szw1710

Kwadratury interpolacyjne

Post autor: szw1710 »

Konstruowanie kwadratur interpolacyjnych
Poniższy tekst powstał w związku z wątkiem https://www.matematyka.pl/viewtopic.php?t=270679 i moją odpowiedzią na jedno z zadanych pytań. Ponieważ jest długi, na pewno wkradły się do niego pewne nieścisłości, błędy w rachunkach, oznaczeniach itp. Bardzo proszę o ich sygnalizowanie na PW oraz o wszelkie uwagi krytyczne.


Przypuśćmy, że dana jest całkowalna funkcja nieujemna \(\displaystyle{ p:[a,b]\to\mathbb{R}}\), dodatnia prawie wszędzie (tzn. \(\displaystyle{ p(x)>0}\) dla prawie wszystkich \(\displaystyle{ x\in[a,b]}\), inaczej, zbiór \(\displaystyle{ \left\{x:p(x)\le 0\right\}}\) ma miarę Lebesgue'a zero). Wtedy \(\displaystyle{ \int_a^bp(x)\text{d}x>0.}\) Funkcję \(\displaystyle{ p}\) nazywamy funkcją wagową.

W praktycznej realizacji wystarczy myśleć, że \(\displaystyle{ p}\) jest funkcją ciągłą, dodatnią, dopuszczając ewentualnie \(\displaystyle{ p(a)=0}\) czy \(\displaystyle{ p(b)=0,}\) choć niekoniecznie.

Naszym zadaniem będzie przybliżenie całki
\(\displaystyle{ \int_a^bf(x)\cdot p(x)\text{d}x}\)
za pomocą kwadratury interpolacyjnej postaci
\(\displaystyle{ \sum_{k=1}^n\lambda_kf(x_k)\,,}\)
gdzie \(\displaystyle{ x_1,\dots,x_n\in[a,b]}\) są ustalonymi punktami, zwanymi węzłami kwadratury, a liczby \(\displaystyle{ \lambda_1,\dots,\lambda_n}\) są ustalonymi współczynnikami. Ponieważ powyższa kwadratura zawiera \(\displaystyle{ n}\) węzłów, nazwiemy ją kwadraturą \(\displaystyle{ n}\)-punktową.

Kwadratury interpolacyjne projektuje się w ten sposób, aby były dokładne na wielomianach określonego rzędu (wielomian jest rzędu \(\displaystyle{ m}\), jeśli jest stopnia co najwyżej \(\displaystyle{ m}\)). Oznacza to, że równość
\(\displaystyle{ \int_a^bw(x)\cdot p(x)\text{d}x=\sum_{k=1}^n\lambda_kw(x_k)}\)
zachodzi dla każdego wielomianu z góry ustalonego rzędu.

Rozpocznijmy od wykazania, że kwadratura \(\displaystyle{ n}\)-punktowa może być dokładna na wielomianach rzędu maksymalnie \(\displaystyle{ 2n-1}\). W tym celu podamy wielomian \(\displaystyle{ w}\) stopnia \(\displaystyle{ 2n}\), dla którego nie zachodzi powyższa równość.

Niech
\(\displaystyle{ p_n(x)=(x-x_1)\cdots(x-x_n)\,,}\)
a \(\displaystyle{ w(x)=p_n^2(x).}\)

Mamy z jednej strony
\(\displaystyle{ \int_a^bw(x)\cdot p(x)\text{d}x>0\,,}\)
co wynika z własności funkcji wagowej i z tego, że \(\displaystyle{ w(x)>0}\) dla prawie wszystkich \(\displaystyle{ x}\), a z drugiej strony
\(\displaystyle{ \sum_{k=1}^n\lambda_kw(x_k)=0\,,}\)
gdyż \(\displaystyle{ x_1,\dots,x_n}\) są pierwiastkami wielomianu \(\displaystyle{ w}\). Nasza kwadratura nie jest więc dokładna na tym wielomianie.

Skonstruujemy teraz kwadraturę w pewnym sensie optymalną, tj. \(\displaystyle{ n}\)-punktową dokładną na wielomianach rzędu \(\displaystyle{ 2n-1}\), bo wyższego rzędu dokładności nie uzyskamy. Niech ma ona postać
\(\displaystyle{ \int_a^b f(x)p(x)\text{d}x\approx \sum_{k=1}^n\lambda_kf(x_k)\,.}\)
Na razie węzły \(\displaystyle{ x_1,\dots,x_n\in[a,b]}\) i współczynniki \(\displaystyle{ \lambda_1,\dots,\lambda_n}\) są bliżej nieokreślone. Spróbujemy dowiedzieć się o nich czegoś więcej.

Najpierw załóżmy, że kwadratura po prawej stronie powyższego wzoru jest dokładna na wielomianach rzędu \(\displaystyle{ 2n-1}\).

W tym miejscu należy wprowadzić parę pojęć.

Rozpocznijmy od całkowego iloczynu skalarnego z funkcją wagową \(\displaystyle{ p}\). Jest to funkcjonał postaci
\(\displaystyle{ \langle f,g\rangle=\int_a^bf(x)g(x)\cdot p(x)\text{d}x}\)
określony na (rzeczywistej) przestrzeni liniowej wszystkich funkcji całkowalnych w przedziale \(\displaystyle{ [a,b].}\) Czytelnikowi pozostawiam sprawdzenie, jest to rzeczywiście iloczyn skalarny. Wystarczy sprawdzić aksjomaty tego iloczynu.

Powiemy, że funkcje całkowalne \(\displaystyle{ f,g:[a,b]\to\mathbb{R}}\)ortogonalne, jeśli \(\displaystyle{ \langle f,g\rangle=0.}\)

Dalej, ciąg wielomianów \(\displaystyle{ (Q_1,Q_2,\dots,Q_n,\dots)}\), gdzie każdy wielomian \(\displaystyle{ Q_k}\) jest stopnia \(\displaystyle{ k}\), jest ortogonalny, jeśli wielomiany \(\displaystyle{ Q_i,Q_j}\) są ortogonalne dla wszystkich \(\displaystyle{ i\ne j.}\) Można zauważyć, że wielomian \(\displaystyle{ Q_k}\) jest wtedy ortogonalny do każdego wielomianu stopnia niższego niż \(\displaystyle{ k}\). Będziemy mówili krótko, że wielomiany \(\displaystyle{ Q_1,Q_2,\dots}\) są ortogonalne w \(\displaystyle{ [a,b]}\) z funkcją wagową \(\displaystyle{ p.}\)

Teoria wielomianów ortogonalnych mówi, że wszystkie wielomiany ortogonalne są wyznaczone jednoznacznie z dokładnością do współczynnika wiodącego (tego przy najwyższej potędze), każdy z nich ma tyle pierwiastków, ile wynosi jego stopień, pierwiastki te są rzeczywiste, jednokrotne i leżą wewnątrz przedziału \(\displaystyle{ [a,b],}\) czyli w przedziale \(\displaystyle{ (a,b).}\)

Weźmy teraz dowolny wielomian \(\displaystyle{ q}\) rzędu \(\displaystyle{ n-1}\). Zatem iloczyn \(\displaystyle{ p_nq}\) jest wielomianem rzędu \(\displaystyle{ 2n-1}\), skąd otrzymujemy
\(\displaystyle{ \langle p_n,q\rangle=\int_a^b p_n(x)q(x)\cdot p(x)\text{d}x=\sum_{k=1}^n\lambda_kp_n(x_k)q(x_k)=0\,,}\)
skąd wynika, że \(\displaystyle{ p_n}\) jest wielomianem ortogonalnym stopnia \(\displaystyle{ n}\) w rodzinie wielomianów ortogonalnych. Wynika stąd, że węzły \(\displaystyle{ x_1,\dots,x_n}\) są pierwiastkami wielomianu ortogonalnego stopnia \(\displaystyle{ n}\).

Wiemy już dużo o węzłach kwadratury \(\displaystyle{ n}\)-punktowej dokładnej na wielomianach rzędu \(\displaystyle{ 2n-1}\). Są one bardzo konkretne. Maksymalny rząd dokładności kwadratury spowodował, że węzły muszą być pierwiastkami wielomianu ortogonalnego.

A co ze współczynnikami \(\displaystyle{ \lambda_1,\dots,\lambda_n}\)?

Żeby je wyznaczyć, konieczne jest następne pojęcie.

Wielomianami podstawowymi Lagrange'a nazywamy wielomiany
\(\displaystyle{ \ell_k(x)=\frac{p_n(x)}{(x-x_k)p_n'(x_k)}\,,\quad k=1,\dots,n\,.}\)
Są one stopnia \(\displaystyle{ n-1}\). Mają bardzo ważną własność:
\(\displaystyle{ \ell_k(x_k)=1,\,\quad\ell_k(x_i)=0\,,\;k\ne i\,,\;k,i=1,\dots,n\,.}\)
Ponadto dla każdej funkcji \(\displaystyle{ f:[a,b]\to\mathbb{R}}\) wielomian \(\displaystyle{ w}\) rzędu \(\displaystyle{ n-1}\) interpolujący tę funkcję w węzłach \(\displaystyle{ x_1,\dots,x_n}\), tj. spełniający warunki interpolacyjne
\(\displaystyle{ w(x_k)=f(x_k)\,,\quad k=1,\dots,n\,,}\)
ma postać
\(\displaystyle{ w(x)=\sum_{k=1}^n\ell_k(x)f(x_k)\,.}\)
Taki wielomian jest wyznaczony jednoznacznie.

Skoro wielomiany podstawowe Lagrange'a są stopnia \(\displaystyle{ n-1}\), to są też rzędu \(\displaystyle{ 2n-1}\). Rozważana przez nas kwadratura jest dokładna na takich wielomianach. Zatem, dla \(\displaystyle{ j=1,\dots,n\,,}\) mamy
\(\displaystyle{ \int_a^b\ell_j(x)\cdot p(x)\text{d}x=\sum_{k=1}^n\lambda_k\ell_j(x_k)=\lambda_k\,.}\)
Dalej, wielomiany \(\displaystyle{ \ell_j^2}\) są stopnia \(\displaystyle{ 2n-2}\), a więc także rzędu \(\displaystyle{ 2n-1}\). Dlatego
\(\displaystyle{ 0<\int_a^b\ell_j^2(x)\cdot p(x)\text{d}x=\sum_{k=1}^n\lambda_k\ell_j^2(x_k)=\lambda_k\,,}\)
co dowodzi dodatniości współczynników \(\displaystyle{ \lambda_1,\dots,\lambda_n.}\)

Podsumowując tę część rozumowania stwierdzamy, że jeśli
\(\displaystyle{ \int_a^bf(x)\cdot p(x)\text{d}x\approx\sum_{k=1}^n\lambda_kf(x_k)}\)
jest \(\displaystyle{ n}\)-punktową kwadraturą dokładną na wielomianach rzędu \(\displaystyle{ 2n-1,}\) to węzły \(\displaystyle{ x_1,\dots,x_n}\) są pierwiastkami wielomianu ortogonalnego \(\displaystyle{ n}\)-tego stopnia w ciągu wielomianów ortogonalnych w \(\displaystyle{ [a,b]}\) z funkcją wagową \(\displaystyle{ p,}\) a współczynniki \(\displaystyle{ \lambda_1,\dots,\lambda_n}\) są dodatnie i wyrażają się wzorami
\(\displaystyle{ \lambda_k=\int_a^b\ell_k(x)\cdot p(x)\text{d}x\,,\quad k=1,\dots,n\,.}\)
Wykażemy, teraz że jeśli \(\displaystyle{ x_1,\dots,x_n,\lambda_1,\dots,\lambda_n}\) są określone jak powyżej, to kwadratura
\(\displaystyle{ \int_a^bf(x)\cdot p(x)\text{d}x\approx\sum_{k=1}^n\lambda_kf(x_k)}\)
jest dokładna na wielomianach rzędu \(\displaystyle{ 2n-1.}\) W tym celu ustalmy dowolny wielomian \(\displaystyle{ w}\) rzędu \(\displaystyle{ 2n-1.}\) Dzieląc go z resztą przez \(\displaystyle{ p_n(x)=(x-x_1)\cdots(x-x_n)}\) otrzymujemy istnienie wielomianów \(\displaystyle{ q}\) (ilorazu) oraz \(\displaystyle{ r}\) (reszty) rzędu \(\displaystyle{ n-1}\) takich, że
\(\displaystyle{ w(x)=p_n(x)q(x)+r(x)\,.}\)
Ale ponieważ węzły kwadratury są pierwiastkami wielomianu ortogonalnego stopnia \(\displaystyle{ n,}\) to \(\displaystyle{ p_n(x)}\) jest, z dokładnością do współczynnika wiodącego, tym właśnie wielomianem. Jest on w szczególności ortogonalny do każdego wielomianu stopnia niższego niż \(\displaystyle{ n.}\) Dlatego
\(\displaystyle{ \int_a^bp_n(x)q(x)\cdot p(x)\text{d}x=0\,.}\)
Zatem
\(\displaystyle{ \int_a^bw(x)\cdot p(x)\text{d}x=\int_a^bp_n(x)q(x)\cdot p(x)\text{d}x+\int_a^br(x)\cdot p(x)\text{d}x=\int_a^br(x)\cdot p(x)\text{d}x\,.}\)
Ale \(\displaystyle{ r,}\) jako wielomian rzędu \(\displaystyle{ n-1,}\) jest równy swojemu wielomianowi interpolacyjnemu Lagrange'a:
\(\displaystyle{ r(x)=\sum_{k=1}^n\ell_k(x)r(x_k)\,.}\)
Dochodzimy więc do
\(\displaystyle{ \begin{aligned}
\int_a^br(x)\cdot p(x)\text{d}x
&=\int_a^b\sum_{k=1}^n\ell_k(x)r(x_k)\cdot p(x)\text{d}x=\\
&=\sum_{k=1}^nr(x_k)\int_a^b\ell_k(x)\cdot p(x)\text{d}x
=\sum_{k=1}^nr(x_k)\lambda_k\,.
\end{aligned}}\)
Biorąc teraz pod uwagę dzielenie wielomianu \(\displaystyle{ w}\) przez \(\displaystyle{ p_n}\) oraz to, że \(\displaystyle{ p_n(x_k)=0}\), \(\displaystyle{ k=1,\dots,n\,,}\) dostajemy
\(\displaystyle{ w(x_k)=p_n(x_k)q(x_k)+r(x_k)=r(x_k)\,,\quad k=1,\dots,n}\)
oraz
\(\displaystyle{ \int_a^br(x)\cdot p(x)\text{d}x=\sum_{k=1}^n\lambda_kr(x_k)=\sum_{k=1}^n\lambda_kw(x_k)}\)
Ostatecznie
\(\displaystyle{ \int_a^bw(x)\cdot p(x)\text{d}x
=\int_a^br(x)\cdot p(x)\text{d}x
=\sum_{k=1}^n\lambda_kw(x_k)\,,}\)
skąd wynika dokładność naszej kwadratury na wielomianach rzędu \(\displaystyle{ 2n-1.}\)

Skonstruowana kwadratura nazywa się kwadraturą Gaussa. Dla przedziału \(\displaystyle{ [-1,1]}\) i funkcji wagowej \(\displaystyle{ p(x)=1}\) otrzymujemy kwadraturę Gaussa-Legendre'a.

Jednopunktowa kwadratura Gaussa-Legendre'a

\(\displaystyle{ \int_{-1}^1f(x)\text{d}x\approx 2f(0)}\)

Dwupunktowa kwadratura Gaussa-Legendre'a

\(\displaystyle{ \int_{-1}^1f(x)\text{d}x\approx f\biggl(-\frac{\sqrt{3}}{3}\biggr)+f\biggl(\frac{\sqrt{3}}{3}\biggr)}\)

Trzypunktowa kwadratura Gaussa-Legendre'a

\(\displaystyle{ \int_{-1}^1f(x)\text{d}x\approx \frac{5}{9}f\biggl(-\frac{\sqrt{15}}{5}\biggr)+\frac{8}{9}f(0)+\frac{5}{9}f\biggl(\frac{\sqrt{15}}{5}\biggr)}\)

Czteropunktowa kwadratura Gaussa-Legendre'a

\(\displaystyle{ \begin{aligned}
\int_{-1}^1f(x)\text{d}x&\approx
\frac{18+\sqrt{{30}}}{36}
\left[
f\left(-\frac{\sqrt{525-70\sqrt{{30}}}}{35}\right)
+f\left(\frac{\sqrt{525-70\sqrt{{30}}}}{35}\right)
\right]+\\[2ex]
&+\frac{18-\sqrt{{30}}}{36}
\left[
f\left(-\frac{\sqrt{525+70\sqrt{{30}}}}{35}\right)
+f\left(\frac{\sqrt{525+70\sqrt{{30}}}}{35}\right)
\right]
\end{aligned}}\)


Pięciopunktowa kwadratura Gaussa-Legendre'a

\(\displaystyle{ \begin{aligned}
\int_{-1}^1f(x)\text{d}x&\approx
\frac{128}{225}f(0)+
\frac{322+13\sqrt{{70}}}{900}
\left[
f\left(-\frac{\sqrt{245-14\sqrt{{70}}}}{21}\right)
+f\left(\frac{\sqrt{245-14\sqrt{{70}}}}{21}\right)
\right]+\\[2ex]
&+\frac{322-13\sqrt{{70}}}{900}
\left[
f\left(-\frac{\sqrt{245+14\sqrt{{70}}}}{21}\right)
+f\left(\frac{\sqrt{245+14\sqrt{{70}}}}{21}\right)
\right]
\end{aligned}}\)


Węzły \(\displaystyle{ x_k}\) kwadratur Gaussa-Legendre'a więcej niż pięciopunktowych (tj. pierwiastki odpowienich wielomianów ortogonalnych zwanych wielomianami Legendre'a) nie wyrażają się przez pierwiastniki, tzn. w formie podobnej do powyższych wzorów. Występuje tu tzw. casus irreducibilis. Podobnie jest ze współczynnikami \(\displaystyle{ \lambda_k}\), gdyż wzory na te współczynniki zawierają węzły \(\displaystyle{ x_k.}\)

Na zakończenie odpowiedź na pytanie dlaczego omawiana kwadratura nazywa się interpolacyjną. Wartość tej kwadratury dla funkcji \(\displaystyle{ f}\), czyli wielkość
\(\displaystyle{ \sum_{k=1}^n\lambda_kf(x_k)\,,}\)
jest równa całce wielomianu \(\displaystyle{ w}\) rzędu \(\displaystyle{ n-1}\) interpolującego funkcję \(\displaystyle{ f}\) w węzłach \(\displaystyle{ x_1,\dots,x_n\,.}\) Można to łatwo sprawdzić. Zapiszmy w tym celu wspomniany wielomian stosując wzór interpolacyjny Lagrange'a:
\(\displaystyle{ w(x)=\sum_{k=1}^n\ell_k(x)f(x_k)}\)
i scałkujmy go:
\(\displaystyle{ \int_a^bw(x)\cdot p(x)\text{d}x=\sum_{k=1}^n\int_a^b\ell_k(x)\cdot p(x)\text{d}x\cdot f(x_k)=\sum_{k=1}^n\lambda_kf(x_k)\,.}\)
Zatem
\(\displaystyle{ \int_a^bf(x)\cdot p(x)\text{d}x
\approx\sum_{k=1}^n\lambda_kf(x_k)=\int_a^bw(x)\cdot p(x)\text{d}x\,.}\)
Można wykazać jeszcze więcej. Tak naprawdę wielkość
\(\displaystyle{ \sum_{k=1}^n\lambda_kf(x_k)\,,}\)
jest równa całce wielomianu \(\displaystyle{ w}\) rzędu \(\displaystyle{ 2n-1}\) interpolującego funkcję \(\displaystyle{ f}\) w węzłach \(\displaystyle{ x_1,\dots,x_n}\) wraz z wartościami pochodnych, tzn. spełniającego warunki interpolacyjne
\(\displaystyle{ w(x_k)=f(x_k)\,,\quad w'(x_k)=f'(x_k)\,,\quad k=1,\dots,n\,.}\)
Wielomian \(\displaystyle{ w}\) zapisujemy stosując wzór interpolacyjny Hermite'a i rozumując analogicznie jak powyżej. Pozwolę sobie pominąć szczegółowe rachunki, które są dość żmudne.
Zablokowany