diff options
author | Andreas Müller <andreas.mueller@ost.ch> | 2021-11-29 13:18:22 +0100 |
---|---|---|
committer | Andreas Müller <andreas.mueller@ost.ch> | 2021-11-29 13:18:22 +0100 |
commit | 02b31342e0dba3703c1a3a91352f7ae19764d7ce (patch) | |
tree | d8463fb987732bd07ef67d2e0a365629af0fe532 | |
parent | gauss quadratur (diff) | |
download | SeminarSpezielleFunktionen-02b31342e0dba3703c1a3a91352f7ae19764d7ce.tar.gz SeminarSpezielleFunktionen-02b31342e0dba3703c1a3a91352f7ae19764d7ce.zip |
gauss quadratur zeugs
-rw-r--r-- | buch/chapters/060-integral/Makefile.inc | 1 | ||||
-rw-r--r-- | buch/chapters/060-integral/gaussquadratur.tex | 484 | ||||
-rw-r--r-- | buch/chapters/060-integral/gq/Makefile | 22 | ||||
-rw-r--r-- | buch/chapters/060-integral/gq/gq.cpp | 64 | ||||
-rw-r--r-- | buch/chapters/060-integral/gq/gq.pdf | bin | 0 -> 25439 bytes | |||
-rw-r--r-- | buch/chapters/060-integral/gq/gq.tex | 57 | ||||
-rw-r--r-- | buch/chapters/060-integral/gq/kreis.cpp | 66 | ||||
-rw-r--r-- | buch/chapters/060-integral/orthogonal.tex | 816 | ||||
-rw-r--r-- | buch/chapters/references.bib | 2 |
9 files changed, 1181 insertions, 331 deletions
diff --git a/buch/chapters/060-integral/Makefile.inc b/buch/chapters/060-integral/Makefile.inc index 1485e40..9cc5356 100644 --- a/buch/chapters/060-integral/Makefile.inc +++ b/buch/chapters/060-integral/Makefile.inc @@ -9,4 +9,5 @@ CHAPTERFILES = $(CHAPTERFILES) \ chapters/060-integral/differentialkoerper.tex \ chapters/060-integral/risch.tex \ chapters/060-integral/orthogonal.tex \ + chapters/060-integral/gaussquadratur.tex \ chapters/060-integral/chapter.tex diff --git a/buch/chapters/060-integral/gaussquadratur.tex b/buch/chapters/060-integral/gaussquadratur.tex new file mode 100644 index 0000000..81d6c28 --- /dev/null +++ b/buch/chapters/060-integral/gaussquadratur.tex @@ -0,0 +1,484 @@ +% +% Anwendung: Gauss-Quadratur +% +\subsection{Anwendung: Gauss-Quadratur} +Orthogonale Polynome haben eine etwas unerwartet Anwendung in einem +von Gauss erdachten numerischen Integrationsverfahren. +Es basiert auf der Beobachtung, dass viele Funktionen sich sehr +gut durch Polynome approximieren lassen. +Wenn man also sicherstellt, dass ein Verfahren für Polynome +sehr gut funktioniert, darf man auch davon ausgehen, dass es für +andere Funktionen nicht allzu schlecht sein wird. + +\subsubsection{Interpolationspolynome} +Zu einer stetigen Funktion $f(x)$ auf dem Intervall $[-1,1]$ +ist ein Polynome vom Grad $n$ gesucht, welches in den Punkten +$x_0<x_1<\dots<x_n$ die Funktionswerte $f(x_i)$ annimmt. +Ein solches Polynom $p(x)$ hat $n+1$ Koeffizienten, die aus dem +linearen Gleichungssystem der $n+1$ Gleichungen $p(x_i)=f(x_i)$ +ermittelt werden können. + +Das Interpolationspolynom $p(x)$ lässt sich abera uch direkt +angeben. +Dazu konstruiert man zuerst die Polynome +\[ +l_i(x) += +\frac{ +(x-x_0)(x-x_1)\cdots\widehat{(x-x_i)}\cdots (x-x_n) +}{ +(x_i-x_0)(x_i-x_1)\cdots\widehat{(x_i-x_i)}\cdots (x_i-x_n) +} +\] +vom Grad $n$, wobei der Hut bedeutet, dass diese Faktoren +im Produkt wegzulassen sind. +Die Polynome $l_i(x)$ haben die Eigenschaft +\[ +l_i(x_j) = \delta_{ij} += +\begin{cases} +1&\qquad i=j\\ +0&\qquad\text{sonst}. +\end{cases} +\] +Die Linearkombination +\[ +p(x) = \sum_{i=0}^n f(x_i)l_i(x) +\] +ist dann ein Polynom vom Grad $n$, welches am den Stellen $x_j$ +die Werte +\[ +p(x_j) += +\sum_{i=0}^n f(x_i)l_i(x_j) += +\sum_{i=0}^n f(x_i)\delta_{ij} += +f(x_j) +\] +hat, das Polynome $p(x)$ ist also das gesuchte Interpolationspolynom. + +\subsubsection{Integrationsverfahren auf der Basis von Interpolation} +Das Integral einer stetigen Funktion $f(x)$ auf dem Intervall $[-1,1]$ +kann mit Hilfe des Interpolationspolynoms approximiert werden. +Wenn $|f(x)-p(x)|<\varepsilon$ ist im Intervall $[-1,1]$, dann gilt +für die Integrale +\[ +\biggl|\int_{-1}^1 f(x)\,dx -\int_{-1}^1p(x)\,dx\biggr| +\le +\int_{-1}^1 |f(x)-p(x)|\,dx +\le +2\varepsilon. +\] +Ein Interpolationspolynom mit kleinem Fehler liefert also auch +eine gute Approximation für das Integral. + +Da das Interpolationspolynome durch die Funktionswerte $f(x_i)$ +bestimmt ist, muss auch das Integral allein aus diesen Funktionswerten +berechnet werden können. +Tatsächlich ist +\begin{equation} +\int_{-1}^1 p(x)\,dx += +\int_{-1}^1 \sum_{i=0}^n f(x_i)l_i(x)\,dx += +\sum_{i=0}^n f(x_i) +\underbrace{\int_{-1}^1 +l_i(x)\,dx}_{\displaystyle = A_i}. +\label{buch:integral:gaussquadratur:eqn:Aidef} +\end{equation} +Das Integral von $f(x)$ wird also durch eine mit den Zahlen $A_i$ +gewichtete Summe +\[ +\int_{-1}^1 f(x)\,dx +\approx +\sum_{i=1}^n f(x_i)A_i +\] +approximiert. + +\subsubsection{Integrationsverfahren, die für Polynome exakt sind} +Ein Polynom vom Grad $2n$ hat $2n+1$ Koeffizienten. +Um das Polynom durch ein Interpolationspolynom exakt wiederzugeben, +braucht man $2n+1$ Stützstellen. +Andererseits gilt +\[ +\int_{-1}^1 a_{2n}x^{2n} + a_{2n-1}x^{2n-1} + \dots + a_2x^2 + a_1x a_0\,dx += +\int_{-1}^1 a_{2n}x^{2n} + a_{2n-2}x^{2n-2}+\dots +a_2x^2 +a_0\,dx, +\] +das Integral ist also bereits durch die $n+1$ Koeffizienten mit geradem +Index bestimmt. +Es sollte daher möglich sein, aus $n+1$ Funktionswerten eines beliebigen +Polynoms vom Grad höchstens $2n$ an geeignet gewählten Stützstellen das +Integral exakt zu bestimmen. + +\begin{beispiel} +Wir versuchen dies für quadratische Polynome durchzuführen, also +für $n=1$. +Gesucht sind also zwei Werte $x_i$, $i=0,1$ und Gewichte $A_i$, $i=0,1$ +derart, dass für jedes quadratische Polynome $p(x)=a_2x^2+a_1x+a_0$ +das Integral durch +\[ +\int_{-1}^1 p(x)\,dx += +A_0 p(x_0) + A_1 p(x_1) +\] +gebeben ist. +Indem wir für $p(x)$ die Polynome $1$, $x$, $x^2$ und $x^3$ einsetzen, +erhalten wir vier Gleichungen +\[ +\begin{aligned} +p(x)&=\rlap{$1$}\phantom{x^2}\colon& 2 &= A_0\phantom{x_0}+ A_1 \\ +p(x)&=x^{\phantom{2}}\colon& 0 &= A_0x_0 + A_1x_1 \\ +p(x)&=x^2\colon& \frac23 &= A_0x_0^2 + A_1x_1^2\\ +p(x)&=x^3\colon& 0 &= A_0x_0^3 + A_1x_1^3. +\end{aligned} +\] +Dividiert man die zweite und vierte Gleichung in der Form +\[ +\left. +\begin{aligned} +A_0x_0 &= -A_1x_1\\ +A_0x_0^2 &= -A_1x_1^2 +\end{aligned} +\quad +\right\} +\quad +\Rightarrow +\quad +x_0^2=x_1^2 +\quad +\Rightarrow +\quad +x_1=-x_0. +\] +Indem wir dies in die zweite Gleichung einsetzen, finden wir +\[ +0 = A_0x_0 + A_1x_1 = A_0x_1 -A_1x_0 = (A_0-A_1)x_0 +\quad\Rightarrow\quad +A_0=A_1. +\] +Aus der ersten Gleichung folgt jetzt +\[ +2= A_0+A_1 = 2A_0 \quad\Rightarrow\quad A_0 = 1. +\] +Damit bleiben nur noch die Werte von $x_i$ zu bestimmen, was +mit Hilfe der zweiten Gleichung geschehen kann: +\[ +\frac23 = A_0x_0^2 + A_1x_1^2 = 2x_0^2 +\quad\Rightarrow\quad +x_0 = \frac{1}{\sqrt{3}}, x_1 = -\frac{1}{\sqrt{3}} +\] +Damit ist das Problem gelöst: das Integral eines Polynoms vom Grad 3 +im Interval $[-1,1]$ ist exakt gegeben durch +\[ +\int_{-1}^1 p(x)\,dx += +p\biggl(-\frac{1}{\sqrt{3}}\biggr) ++ +p\biggl(\frac{1}{\sqrt{3}}\biggr). +\] +Das Integral kann also durch nur zwei Auswertungen des Polynoms +exakt bestimmt werden. + +Im Laufe der Lösung des Gleichungssystems wurden die Gewichte $A_i$ +mit bestimmt. +Es ist aber auch möglich, die Gewichte zu bestimmen, wenn man die +Stützstellen kennt. +Nach \eqref{buch:integral:gaussquadratur:eqn:Aidef} +sind sie die $A_i$ gegeben als Integrale der Polynome +$l_i(x)$, die im vorliegenden Fall linear sind: +\begin{align*} +l_0(x) +&= +\frac{x-x_1}{x_0-x_1} += +\frac{x-\frac1{\sqrt{3}}}{-\frac{2}{\sqrt{3}}} += +\frac12(1-\sqrt{3}x) +\\ +l_1(x) +&= +\frac{x-x_0}{x_1-x_0} += +\frac{x+\frac1{\sqrt{3}}}{\frac{2}{\sqrt{3}}} += +\frac12(1+\sqrt{3}x) +\end{align*} +Diese haben die Integrale +\[ +\int_{-1}^1\frac12(1\pm\sqrt{3}x)\,dx += +\int_{-1}^1 \frac12\,dx += +1, +\] +da das Polynom $x$ verschwindendes Integral hat. +Dies stimmt mit $A_0=A_1=1$ überein. +\label{buch:integral:beispiel:gaussquadraturn1} +\end{beispiel} + +Das eben vorgestellt Verfahren kann natürlich auf beliebiges $n$ +verallgemeinert werden. +Allerdings ist die Rechnung zur Bestimmung der Stützstellen und +Gewichte sehr mühsam. + +\subsubsection{Stützstellen und Orthogonalpolynome} +Sei $R_n=\{p(X)\in\mathbb{R}[X] \mid \deg p\le n\}$ der Vektorraum +der Polynome vom Grad $n$. + +\begin{satz} +\label{buch:integral:satz:gaussquadratur} +Sei $p$ ein Polynom vom Grad $n$, welches auf allen Polynomen in $R_{n-1}$ +orthogonal sind. +Seien ausserdem $x_0<x_1<\dots<x_n$ Stützstellen im Intervall $[-1,1]$ +und $A_i\in\mathbb{R}$ Gewichte derart dass +\[ +\int_{-1}^1 f(x)\,dx = +\sum_{i=0}^n A_if(x_i) +\] +für jedes Polynom $f$ vom Grad höchstens $2n-1$, dann sind die Zahlen +$x_i$ die Nullstellen des Polynoms $p$. +\end{satz} + +\begin{proof}[Beweis] +Sei $f(x)$ ein beliebiges Polynom vom Grad $2n-1$. +Nach dem Polynomdivisionsalgorithmus gibt es +Polynome $q,r\in R_{n-1}$ derart, dass $f=qp+r$. +Dann ist das Integral von $f$ gegeben durch +\[ +\int_{-1}^1 f(x)\,dx += +\int_{-1}^1q(x) p(x)\,dx + \int_{-1}^1 r(x)\,dx += +\langle q,p\rangle + \int_{-1}^1 r(x)\,dx. +\] +Da $p\perp R_{n-1}$ folgt insbesondere, dass $\langle q,p\rangle=0$. + +Da die Integrale auch aus den Werten in den Stützstellen berechnet +werden können, muss auch +\[ +0 += +\int_{-1}^1 q(x)p(x)\,dx += +\sum_{i=0}^n q(x_i)p(x_i) +\] +für jedes beliebige Polynom $q\in R_{n-1}$ gelten. +Da man für $q$ die Interpolationspolynome $l_j(x)$ verwenden +kann, den Grad $n-1$ haben, folgt +\[ +0 += +\sum_{i=0}^n +l_j(x_i)p(x_i) += +\sum_{i=0}^n \delta_{ij}p(x_i), +\] +die Stützstellen $x_i$ müssen also die Nullstellen des Polynoms +$p(x)$ sein. +\end{proof} + +Der Satz~\ref{buch:integral:satz:gaussquadratur} begründet das +{\em Gausssche Quadraturverfahren}. +Die in Abschnitt~\ref{buch:integral:subsection:orthogonale-polynome} +bestimmten Legendre-Polynome $P_n$ haben die im Satz +verlangte Eigenschaft, +dass sie auf allen Polynomen geringeren Grades orthogonal sind. +Wählt man die $n$ Nullstellen von $P_n$ als Stützstellen, erhält man +automatisch ein Integrationsverfahren, welches für Polynome vom Grad +$2n-1$ exakt ist. + +\begin{beispiel} +Das Legendre-Polynom $P_2(x) = \frac12(3x^2-1)$ hat die +Nullstellen $x=\pm1/\sqrt{3}$, dies sind genau die im Beispiel +auf Seite~\pageref{buch:integral:beispiel:gaussquadraturn1} befundenen +Sützstellen. +\end{beispiel} + +\subsubsection{Fehler der Gauss-Quadratur} +Das Gausssche Quadraturverfahren mit $n$ Stützstellen berechnet +Integrale von Polynomen bis zum Grad $2n-1$ exakt. +Für eine beliebige Funktion kann man die folgende Fehlerabschätzung +angeben \cite[theorem 7.3.4, p.~497]{buch:numal}. + +\begin{satz} +Seien $x_i$ die Stützstellen und $A_i$ die Gewichte einer +Gaussschen Quadraturformel mit $n+1$ Stützstellen und sei $f$ +eine auf dem Interval $[-1,1]$ $2n+2$-mal stetig differenzierbare +Funktion, dann ist der $E$ Fehler des Integrals +\[ +\int_{-1}^1 f(x)\,dx = \sum_{i=0}^n A_i f(x_i) + E +\] +gegeben durch +\begin{equation} +E = \frac{f^{(2n+2)}(\xi)}{(2n+2)!}\int_{-1}^1 l(x)^2\,dx, +\label{buch:integral:gaussquadratur:eqn:fehlerformel} +\end{equation} +wobei $l(x)=(x-x_0)(x-x_1)\dots(x-x_n)$ und $\xi$ ein geeigneter +Wert im Intervall $[-1,1]$ ist. +\end{satz} + +Dank dem Faktor $(2n+2)!$ im Nenner von +\eqref{buch:integral:gaussquadratur:eqn:fehlerformel} +geht der Fehler für grosses $n$ sehr schnell gegen $0$. +Man kann auch zeigen, dass die mit Gauss-Quadratur mit $n+1$ +Stützstellen berechneten Näherungswerte eines Integrals einer +stetigen Funktion $f(x)$ für $n\to\infty$ immer gegen den wahren +Wert des Integrals konvergieren. + +\begin{table} +\def\u#1{\underline{#1}} +\centering +\begin{tabular}{|>{$}c<{$}|>{$}r<{$}|>{$}r<{$}|} +\hline + n & \text{Gauss-Quadratur} & \text{Trapezregel} \\ +\hline +\phantom{0}2 & 0.\u{95}74271077563381 & 0.\u{95}63709682242596 \\ +\phantom{0}4 & 0.\u{95661}28333449730 & 0.\u{956}5513401768598 \\ +\phantom{0}6 & 0.\u{9566114}812034364 & 0.\u{956}5847489712136 \\ +\phantom{0}8 & 0.\u{956611477}5028123 & 0.\u{956}5964425360520 \\ + 10 & 0.\u{9566114774905}637 & 0.\u{9566}018550715587 \\ + 12 & 0.\u{956611477490518}7 & 0.\u{9566}047952369826 \\ + 14 & 0.\u{95661147749051}72 & 0.\u{9566}065680717177 \\ + 16 & 0.\u{956611477490518}7 & 0.\u{9566}077187127541 \\ + 18 & 0.\u{956611477490518}3 & 0.\u{9566}085075898731 \\ + 20 & 0.\u{956611477490518}4 & 0.\u{9566}090718697414 \\ +\hline + \infty & 0.9566114774905183 & 0.9566114774905183 \\ +\hline +\end{tabular} +\caption{Integral von $\sqrt{1-x^2}$ zwischen $-\frac12$ und $\frac12$ +berechnet mit Gauss-Quadratur und der Trapezregel, aber mit zehnmal +so vielen Stützstellen. +Bereits mit 12 Stützstellen erreicht die Gauss-Quadratur +Maschinengenauigkeit, die Trapezregel liefert auch mit 200 Stützstellen +nicht mehr als 4 korrekte Nachkommastellen. +\label{buch:integral:gaussquadratur:table0.5}} +\end{table} + +%\begin{table} +%\def\u#1{\underline{#1}} +%\centering +%\begin{tabular}{|>{$}c<{$}|>{$}r<{$}|>{$}r<{$}|} +%\hline +% n & \text{Gauss-Quadratur} & \text{Trapezregel} \\ +%\hline +%\phantom{0}2 & 1.\u{5}379206741571556 & 1.\u{5}093105464758343 \\ +%\phantom{0}4 & 1.\u{51}32373472933831 & 1.\u{51}13754509594814 \\ +%\phantom{0}6 & 1.\u{512}1624557410367 & 1.\u{51}17610879524799 \\ +%\phantom{0}8 & 1.\u{51207}93479994321 & 1.\u{51}18963282632112 \\ +% 10 & 1.\u{51207}13859966004 & 1.\u{51}19589735776959 \\ +% 12 & 1.\u{512070}5317779943 & 1.\u{51}19930161260693 \\ +% 14 & 1.\u{5120704}334802813 & 1.\u{5120}135471596636 \\ +% 16 & 1.\u{5120704}216176006 & 1.\u{5120}268743889558 \\ +% 18 & 1.\u{5120704}201359081 & 1.\u{5120}360123137213 \\ +% 20 & 1.\u{5120704199}459651 & 1.\u{5120}425490275837 \\ +%\hline +% \infty & 1.5120704199172947 & 1.5120704199172947 \\ +%\hline +%\end{tabular} +%\end{table} + +%\begin{table} +%\def\u#1{\underline{#1}} +%\centering +%\begin{tabular}{|>{$}c<{$}|>{$}r<{$}|>{$}r<{$}|} +%\hline +% n & \text{Gauss-Quadratur} & \text{Trapezregel} \\ +%\hline +%\phantom{0}2 & 1.\u{}6246862220133462 & 1.\u{5}597986803933712 \\ +%\phantom{0}4 & 1.\u{5}759105515463101 & 1.\u{56}63563456168101 \\ +%\phantom{0}6 & 1.\u{5}706630058381434 & 1.\u{56}77252866190838 \\ +%\phantom{0}8 & 1.\u{56}94851106536780 & 1.\u{568}2298707696152 \\ +% 10 & 1.\u{56}91283195332679 & 1.\u{568}4701957758742 \\ +% 12 & 1.\u{56}90013806299465 & 1.\u{568}6030805941198 \\ +% 14 & 1.\u{5689}515434853885 & 1.\u{568}6841603070025 \\ +% 16 & 1.\u{5689}306507843050 & 1.\u{568}7372230731711 \\ +% 18 & 1.\u{5689}214761291217 & 1.\u{568}7738235496322 \\ +% 20 & 1.\u{56891}73062385982 & 1.\u{568}8001228530786 \\ +%\hline +% \infty & 1.5689135396691616 & 1.5689135396691616 \\ +%\hline +%\end{tabular} +%\end{table} + +\begin{table} +\def\u#1{\underline{#1}} +\centering +\begin{tabular}{|>{$}c<{$}|>{$}r<{$}|>{$}r<{$}|} +\hline + n & \text{Gauss-Quadratur} & \text{Trapezregel} \\ +\hline +\phantom{0}2 & 1.\u{}6321752373234928 & 1.\u{5}561048774629949 \\ +\phantom{0}4 & 1.\u{57}98691557134743 & 1.\u{5}660124134617943 \\ +\phantom{0}6 & 1.\u{57}35853681692993 & 1.\u{5}683353001877542 \\ +\phantom{0}8 & 1.\u{57}19413565928206 & 1.\u{5}692627503425400 \\ + 10 & 1.\u{57}13388119633434 & 1.\u{5}697323578543481 \\ + 12 & 1.\u{57}10710489948883 & 1.\u{570}0051217458713 \\ + 14 & 1.\u{570}9362135398341 & 1.\u{570}1784766276063 \\ + 16 & 1.\u{570}8621102742815 & 1.\u{570}2959121005231 \\ + 18 & 1.\u{570}8186779483588 & 1.\u{570}3793521168343 \\ + 20 & 1.\u{5707}919411931615 & 1.\u{570}4408749735932 \\ +\hline + \infty & 1.5707367072605671 & 1.5707367072605671 \\ +\hline +\end{tabular} +\caption{Integral von $\sqrt{1-x^2}$ zwischen $-0.999$ und $0.999$ +berechnet mit Gauss-Quadratur und der Trapezregel, aber mit zehnmal +so vielen Stützstellen. +Wegen der divergierenden Steigung des Integranden bei $\pm 1$ tun +sich beide Verfahren sehr schwer. +Trotzdem erreich die Gauss-Quadrator 4 korrekte Nachkommastellen +mit 20 Stütztstellen, während die Trapezregel auch mit 200 Stützstellen +nur 3 korrekte Nachkommastellen findet. +\label{buch:integral:gaussquadratur:table0.999}} +\end{table} + +\begin{figure} +\centering +\includegraphics{chapters/060-integral/gq/gq.pdf} +\caption{Approximationsfehler des +Integrals~\eqref{buch:integral:gaussquadratur:bspintegral} +in Abhängigkeit von $a$. +Die Divergenz der Ableitung des Integranden an den Intervallenden +$\pm 1$ führt zu schlechter Konvergenz des Verfahrens, wenn $a$ +nahe an $1$ ist. +\label{buch:integral:gaussquadratur:fehler}} +\end{figure} + +Zur Illustration der Genauigkeit der Gauss-Quadratur berechnen wir +das Integral +\begin{equation} +\int_{-a}^a \sqrt{1-x^2}\,dx += +\arcsin a + a \sqrt{1-a^2} +\label{buch:integral:gaussquadratur:bspintegral} +\end{equation} +mit Gauss-Quadratur einerseits und dem Trapezverfahren +andererseits. +Da Gauss-Quadratur mit sehr viel weniger Sützstellen auskommt, +berechnen wir die Trapeznäherung mit zehnmal so vielen Stützstelln. +In den Tabellen~\ref{buch:integral:gaussquadratur:table0.5} +und +\ref{buch:integral:gaussquadratur:table0.999} +sind die Resultate zusammengestellt. +Für $a =\frac12$ zeigt +Tabelle~\ref{buch:integral:gaussquadratur:table0.5} +die sehr schnelle Konvergenz der Gauss-Quadratur, schon mit +12 Stützstellen wird Maschinengenauigkeit erreicht. +Das Trapezverfahren dagegen erreicht auch mit 200 Stützstellen nur +4 korrekte Nachkommastellen. + +An den Stellen $x=\pm 1$ divergiert die Ableitung des Integranden +des Integrals \eqref{buch:integral:gaussquadratur:bspintegral}. +Da grösste und kleinste Stützstelle der Gauss-Quadratur immer +deutlich vom Rand des Intervalls entfernt ist, kann das Verfahren +diese ``schwierigen'' Stellen nicht erkennen. +Tabelle~\ref{buch:integral:gaussquadratur:table0.999} zeigt, wie +die Konvergenz des Verfahrens in diesem Fall sehr viel schlechter ist. +Dies zeigt auch der Graph in +Abbildung~\ref{buch:integral:gaussquadratur:fehler}. + +\subsubsection{Skalarprodukte mit Gewichtsfunktion} + diff --git a/buch/chapters/060-integral/gq/Makefile b/buch/chapters/060-integral/gq/Makefile new file mode 100644 index 0000000..2ce859c --- /dev/null +++ b/buch/chapters/060-integral/gq/Makefile @@ -0,0 +1,22 @@ +# +# Makefile -- Gauss quadrature experiments +# +# (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule +# +all: kreis gq.pdf + +gq.o: gq.cpp + g++ `pkg-config --cflags gsl` -o gq.o -g -Wall -O2 -c gq.cpp + +gq: gq.o + g++ -g -Wall -O2 gq.o `pkg-config --libs gsl` -o gq + +kreis: kreis.cpp + g++ `pkg-config --cflags gsl` -g -Wall -O2 kreis.cpp \ + `pkg-config --libs gsl` -o kreis + +gqpaths.tex: gq + ./gq > gqpaths.tex + +gq.pdf: gq.tex gqpaths.tex + pdflatex gq.tex diff --git a/buch/chapters/060-integral/gq/gq.cpp b/buch/chapters/060-integral/gq/gq.cpp new file mode 100644 index 0000000..3eb3ac8 --- /dev/null +++ b/buch/chapters/060-integral/gq/gq.cpp @@ -0,0 +1,64 @@ +/* + * gq.cpp -- Gauss Quadrature Example + * + * (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule + */ +#include <iostream> +#include <iomanip> +#include <cstdlib> +#include <cstdio> +#include <cmath> +#include <gsl/gsl_integration.h> + +int counter; + +double f(double x, void *) { + counter++; + return sqrt(1-x*x); +} + +void do_integration(size_t n, const char *name) { + double scaling = 1e6; + // set up integration gauss quadrature + counter = 0; + gsl_function func; + func.function = f; + func.params = NULL; + + // prepare plotting of results + int steps = 200; + double h = 1. / steps; + + printf("\\def\\%s{\n\t(0,0)", name); + for (int i = 1; i <= steps; i++) { + double a = i * h; + a = 1 - (1 - a) * (1 - a); + double expected = (asin(a) + a * sqrt(1 - a*a)) / 2; + gsl_integration_fixed_workspace *ws + = gsl_integration_fixed_alloc( + gsl_integration_fixed_legendre, + n, 0, a, 0, 0); + double result = 0; + gsl_integration_fixed(&func, &result, ws); + double y = fabs(scaling * (result - expected)); + if (y > 10) { + y = 10; + } + printf("\n\t--({\\dx*%.4f},{\\dy*%.6f})", a, y); + gsl_integration_fixed_free(ws); + } + printf("\n}\n"); +} + +int main(int argc, char *argv[]) { + do_integration(2, "two"); + do_integration(4, "four"); + do_integration(6, "six"); + do_integration(8, "eight"); + do_integration(10, "ten"); + do_integration(12, "twelve"); + do_integration(14, "fourteen"); + do_integration(16, "sixteen"); + + return EXIT_SUCCESS; +} diff --git a/buch/chapters/060-integral/gq/gq.pdf b/buch/chapters/060-integral/gq/gq.pdf Binary files differnew file mode 100644 index 0000000..7a0b3c4 --- /dev/null +++ b/buch/chapters/060-integral/gq/gq.pdf diff --git a/buch/chapters/060-integral/gq/gq.tex b/buch/chapters/060-integral/gq/gq.tex new file mode 100644 index 0000000..e5a9c4d --- /dev/null +++ b/buch/chapters/060-integral/gq/gq.tex @@ -0,0 +1,57 @@ +% +% gq.tex -- Fehler +% +% (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule +% +\documentclass[tikz]{standalone} +\usepackage{amsmath} +\usepackage{times} +\usepackage{txfonts} +\usepackage{pgfplots} +\usepackage{csvsimple} +\usetikzlibrary{arrows,intersections,math} +\begin{document} +\def\skala{1} +\def\dx{12} +\def\dy{8} +\input{gqpaths.tex} +\begin{tikzpicture}[>=latex,thick,scale=\skala] + +\begin{scope} +\clip (0,0) rectangle ({\dx*1},{\dy*1}); +\draw[color=red,line width=1.4pt] \two; +\draw[color=red,line width=1.4pt] \four; +\draw[color=red,line width=1.4pt] \six; +\draw[color=red,line width=1.4pt] \eight; +\draw[color=red,line width=1.4pt] \ten; +\draw[color=red,line width=1.4pt] \twelve; +\draw[color=red,line width=1.4pt] \fourteen; +\draw[color=red,line width=1.4pt] \sixteen; +\end{scope} + +% add image content here +\draw[->] (0,-0.1) -- (0,{\dy+0.5}) coordinate[label={right:$e$}]; +\draw[->] (-0.1,0) -- ({\dx+0.5},0) coordinate[label={$a$}]; + +\draw ({\dx},-0.1) -- ({\dx},0.1); +\draw[line width=0.3pt] ({\dx},0) -- ({\dx},{\dy}); +\node at ({\dx},0) [below] {$1$}; +\node at (0,0) [below] {$0$}; + +\draw (-0.1,\dy) -- (0.1,\dy); +\node at (0,\dy) [left] {$10^{-6}$}; + +\node at ({0.227*\dx},{0.7*\dy}) [rotate=83.5] {$n=2$}; +\node at ({0.678*\dx},{0.7*\dy}) [rotate=83] {$n=4$}; +\node at ({0.850*\dx},{0.7*\dy}) [rotate=84] {$n=6$}; +\node at ({0.915*\dx},{0.7*\dy}) [rotate=85] {$n=8$}; +\node at ({1.015*\dx},{0.7*\dy}) [rotate=90] {$n=16$}; + +\foreach \x in {0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9}{ + \draw ({\x*\dx},-0.1) -- ({\x*\dx},0.1); + \node at ({\x*\dx},0) [below] {$\x$}; +} + +\end{tikzpicture} +\end{document} + diff --git a/buch/chapters/060-integral/gq/kreis.cpp b/buch/chapters/060-integral/gq/kreis.cpp new file mode 100644 index 0000000..c9723cc --- /dev/null +++ b/buch/chapters/060-integral/gq/kreis.cpp @@ -0,0 +1,66 @@ +/* + * kreis.cpp -- Gauss Quadrature Example + * + * (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule + */ +#include <iostream> +#include <iomanip> +#include <cstdlib> +#include <cstdio> +#include <cmath> +#include <gsl/gsl_integration.h> + +double xmax = 0.990; + +double f(double x, void *) { + return sqrt(1-x*x); +} + +double do_integration(size_t n) { + // set up integration gauss quadrature + gsl_function func; + func.function = f; + func.params = NULL; + + gsl_integration_fixed_workspace *ws + = gsl_integration_fixed_alloc( + gsl_integration_fixed_legendre, + n, -xmax, xmax, 0, 0); + double result = 0; + gsl_integration_fixed(&func, &result, ws); + gsl_integration_fixed_free(ws); + return result; +} + +double do_trapez(size_t n) { + double h = 2 * xmax / n; + double s = 0; + for (int i = 0; i <= n; i++) { + double x = -xmax + i * h; + double a = f(x, NULL); + if ((i == 0) || (i == n)) { + a = a * 0.5; + } + s += a; + } + return h * s; +} + +void do_table(double limit) { + xmax = limit; + for (int n = 2; n <= 20; n += 2) { + printf("%d & %.16f & %.16f \\\\\n", n, do_integration(n), + do_trapez(10 * n)); + } + double amax = acos(xmax); + double expected = (M_PI / 2 - amax) + sin(amax) * cos(amax); + printf("\\infty & %.16f & %.16f \\\\\n", expected, expected); +} + +int main(int argc, char *argv[]) { + do_table(0.5); + do_table(0.9); + do_table(0.99); + do_table(0.999); + return EXIT_SUCCESS; +} diff --git a/buch/chapters/060-integral/orthogonal.tex b/buch/chapters/060-integral/orthogonal.tex index 4ffbbde..ceba53a 100644 --- a/buch/chapters/060-integral/orthogonal.tex +++ b/buch/chapters/060-integral/orthogonal.tex @@ -5,6 +5,7 @@ % \section{Orthogonalität \label{buch:integral:section:orthogonale-polynome}} +\rhead{Orthogonale Polynome} Die Fourier-Theorie basiert auf der Idee, Funktionen durch Funktionenreihen mit Summanden zu bilden, die im Sinne eines Skalarproduktes orthogonal sind, welches mit Hilfe eines Integrals @@ -294,6 +295,7 @@ sind sie jedoch nicht orthogonal, denn es ist \biggl[\frac{x^{i+j+1}}{i+j+1}\biggr]_{-1}^1 = \begin{cases} +\displaystyle \frac{2}{i+j+1}&\qquad\text{$i+j$ gerade}\\ 0&\qquad\text{$i+j$ ungerade}. \end{cases} @@ -454,7 +456,8 @@ x^4 &= x^4 -\frac47 P_2(x) - \frac15 P_0(x) -= +\\ +&= x^4 - \frac{6}{7}x^2 + \frac{3}{35} \end{align*} mit $p(1)=\frac{8}{35}$, so dass man @@ -514,334 +517,487 @@ Tabelle~\ref{buch:integral:table:legendre-polynome}. \subsubsection{Legendre-Funktionen zweiter Art} Siehe Wikipedia-Artikel \url{https://de.wikipedia.org/wiki/Legendre-Polynom} +%% +%% Anwendung: Gauss-Quadratur +%% +%\subsection{Anwendung: Gauss-Quadratur} +%Orthogonale Polynome haben eine etwas unerwartet Anwendung in einem +%von Gauss erdachten numerischen Integrationsverfahren. +%Es basiert auf der Beobachtung, dass viele Funktionen sich sehr +%gut durch Polynome approximieren lassen. +%Wenn man also sicherstellt, dass ein Verfahren für Polynome +%sehr gut funktioniert, darf man auch davon ausgehen, dass es für +%andere Funktionen nicht allzu schlecht sein wird. % -% Anwendung: Gauss-Quadratur +%\subsubsection{Interpolationspolynome} +%Zu einer stetigen Funktion $f(x)$ auf dem Intervall $[-1,1]$ +%ist ein Polynome vom Grad $n$ gesucht, welches in den Punkten +%$x_0<x_1<\dots<x_n$ die Funktionswerte $f(x_i)$ annimmt. +%Ein solches Polynom $p(x)$ hat $n+1$ Koeffizienten, die aus dem +%linearen Gleichungssystem der $n+1$ Gleichungen $p(x_i)=f(x_i)$ +%ermittelt werden können. % -\subsection{Anwendung: Gauss-Quadratur} -Orthogonale Polynome haben eine etwas unerwartet Anwendung in einem -von Gauss erdachten numerischen Integrationsverfahren. -Es basiert auf der Beobachtung, dass viele Funktionen sich sehr -gut durch Polynome approximieren lassen. -Wenn man also sicherstellt, dass ein Verfahren für Polynome -sehr gut funktioniert, darf man auch davon ausgehen, dass es für -andere Funktionen nicht allzu schlecht sein wird. - -\subsubsection{Interpolationspolynome} -Zu einer stetigen Funktion $f(x)$ auf dem Intervall $[-1,1]$ -ist ein Polynome vom Grad $n$ gesucht, welches in den Punkten -$x_0<x_1<\dots<x_n$ die Funktionswerte $f(x_i)$ annimmt. -Ein solches Polynom $p(x)$ hat $n+1$ Koeffizienten, die aus dem -linearen Gleichungssystem der $n+1$ Gleichungen $p(x_i)=f(x_i)$ -ermittelt werden können. - -Das Interpolationspolynom $p(x)$ lässt sich abera uch direkt -angeben. -Dazu konstruiert man zuerst die Polynome -\[ -l_i(x) -= -\frac{ -(x-x_0)(x-x_1)\cdots\widehat{(x-x_i)}\cdots (x-x_n) -}{ -(x_i-x_0)(x_i-x_1)\cdots\widehat{(x_i-x_i)}\cdots (x_i-x_n) -} -\] -vom Grad $n$, wobei der Hut bedeutet, dass diese Faktoren -im Produkt wegzulassen sind. -Die Polynome $l_i(x)$ haben die Eigenschaft -\[ -l_i(x_j) = \delta_{ij} -= -\begin{cases} -1&\qquad i=j\\ -0&\qquad\text{sonst}. -\end{cases} -\] -Die Linearkombination -\[ -p(x) = \sum_{i=0}^n f(x_i)l_i(x) -\] -ist dann ein Polynom vom Grad $n$, welches am den Stellen $x_j$ -die Werte -\[ -p(x_j) -= -\sum_{i=0}^n f(x_i)l_i(x_j) -= -\sum_{i=0}^n f(x_i)\delta_{ij} -= -f(x_j) -\] -hat, das Polynome $p(x)$ ist also das gesuchte Interpolationspolynom. - -\subsubsection{Fehler des Interpolationspolynoms} -TODO - -\subsubsection{Integrationsverfahren auf der Basis von Interpolation} -Das Integral einer stetigen Funktion $f(x)$ auf dem Intervall $[-1,1]$ -kann mit Hilfe des Interpolationspolynoms approximiert werden. -Wenn $|f(x)-p(x)|<\varepsilon$ ist im Intervall $[-1,1]$, dann gilt -für die Integrale -\[ -\biggl|\int_{-1}^1 f(x)\,dx -\int_{-1}^1p(x)\,dx\biggr| -\le -\int_{-1}^1 |f(x)-p(x)|\,dx -\le -2\varepsilon. -\] -Ein Interpolationspolynom mit kleinem Fehler liefert also auch -eine gute Approximation für das Integral. - -Da das Interpolationspolynome durch die Funktionswerte $f(x_i)$ -bestimmt ist, muss auch das Integral allein aus diesen Funktionswerten -berechnet werden können. -Tatsächlich ist -\[ -\int_{-1}^1 p(x)\,dx -= -\int_{-1}^1 \sum_{i=0}^n f(x_i)l_i(x)\,dx -= -\sum_{i=0}^n f(x_i) -\underbrace{\int_{-1}^1 -l_i(x)\,dx}_{\displaystyle = A_i} -\] -Das Integral von $f(x)$ wird also durch eine mit den Zahlen $A_i$ -gewichtete Summe -\[ -\int_{-1}^1 f(x)\,dx -\approx -\sum_{i=1}^n f(x_i)A_i -\] -approximiert. - -\subsubsection{Integrationsverfahren, die für Polynome exakt sind} -Ein Polynom vom Grad $2n$ hat $2n+1$ Koeffizienten. -Um das Polynom durch ein Interpolationspolynom exakt wiederzugeben, -braucht man $2n+1$ Stützstellen. -Andererseits gilt -\[ -\int_{-1}^1 a_{2n}x^{2n} + a_{2n-1}x^{2n-1} + \dots + a_2x^2 + a_1x a_0\,dx -= -\int_{-1}^1 a_{2n}x^{2n} + a_{2n-2}x^{2n-2}+\dots +a_2x^2 +a_0\,dx, -\] -das Integral ist also bereits durch die $n+1$ Koeffizienten mit geradem -Index bestimmt. -Es sollte daher möglich sein, aus $n+1$ Funktionswerten eines beliebigen -Polynoms vom Grad höchstens $2n$ an geeignet gewählten Stützstellen das -Integral exakt zu bestimmen. - -\begin{beispiel} -Wir versuchen dies für quadratische Polynome durchzuführen, also -für $n=1$. -Gesucht sind also zwei Werte $x_i$, $i=0,1$ und Gewichte $A_i$, $i=0,1$ -derart, dass für jedes quadratische Polynome $p(x)=a_2x^2+a_1x+a_0$ -das Integral durch -\[ -\int_{-1}^1 p(x)\,dx -= -A_0 p(x_0) + A_1 p(x_1) -\] -gebeben ist. -Indem wir für $p(x)$ die Polynome $1$, $x$, $x^2$ und $x^3$ einsetzen, -erhalten wir vier Gleichungen -\[ -\begin{aligned} -p(x)&=\rlap{$1$}\phantom{x^2}\colon& 2 &= A_0\phantom{x_0}+ A_1 \\ -p(x)&=x^{\phantom{2}}\colon& 0 &= A_0x_0 + A_1x_1 \\ -p(x)&=x^2\colon& \frac23 &= A_0x_0^2 + A_1x_1^2\\ -p(x)&=x^3\colon& 0 &= A_0x_0^3 + A_1x_1^3. -\end{aligned} -\] -Dividiert man die zweite und vierte Gleichung in der Form -\[ -\left. -\begin{aligned} -A_0x_0 &= -A_1x_1\\ -A_0x_0^2 &= -A_1x_1^2 -\end{aligned} -\quad -\right\} -\quad -\Rightarrow -\quad -x_0^2=x_1^2 -\quad -\Rightarrow -\quad -x_1=-x_0. -\] -Indem wir dies in die zweite Gleichung einsetzen, finden wir -\[ -0 = A_0x_0 + A_1x_1 = A_0x_1 -A_1x_0 = (A_0-A_1)x_0 -\quad\Rightarrow\quad -A_0=A_1. -\] -Aus der ersten Gleichung folgt jetzt -\[ -2= A_0+A_1 = 2A_0 \quad\Rightarrow\quad A_0 = 1. -\] -Damit bleiben nur noch die Werte von $x_i$ zu bestimmen, was -mit Hilfe der zweiten Gleichung geschehen kann: -\[ -\frac23 = A_0x_0^2 + A_1x_1^2 = 2x_0^2 -\quad\Rightarrow\quad -x_0 = \frac{1}{\sqrt{3}}, x_1 = -\frac{1}{\sqrt{3}} -\] -Damit ist das Problem gelöst: das Integral eines Polynoms vom Grad 3 -im Interval $[-1,1]$ ist exakt gegeben durch -\[ -\int_{-1}^1 p(x)\,dx -= -p\biggl(-\frac{1}{\sqrt{3}}\biggr) -+ -p\biggl(\frac{1}{\sqrt{3}}\biggr). -\] -Das Integral kann also durch nur zwei Auswertungen des Polynoms -exakt bestimmt werden. - -Im Laufe der Lösung des Gleichungssystems wurden die Gewichte $A_i$ -mit bestimmt. -Es ist aber auch möglich, die Gewichte zu bestimmen, wenn man die -Stützstellen kennt. -Nach \eqref{XXX} sind sie gegeben als Integrale der Polynome $l_i(x)$, -die im vorliegenden Fall linear sind: -\begin{align*} -l_0(x) -&= -\frac{x-x_1}{x_0-x_1} -= -\frac{x-\frac1{\sqrt{3}}}{-\frac{2}{\sqrt{3}}} -= -\frac12(1-\sqrt{3}x) -\\ -l_1(x) -&= -\frac{x-x_0}{x_1-x_0} -= -\frac{x+\frac1{\sqrt{3}}}{\frac{2}{\sqrt{3}}} -= -\frac12(1+\sqrt{3}x) -\end{align*} -Diese haben die Integrale -\[ -\int_{-1}^1\frac12(1\pm\sqrt{3}x)\,dx -= -\int_{-1}^1 \frac12\,dx -= -1, -\] -da das Polynom $x$ verschwindendes Integral hat. -Dies stimmt mit $A_0=A_1=1$ überein. -\label{buch:integral:beispiel:gaussquadraturn1} -\end{beispiel} - -Das eben vorgestellt Verfahren kann natürlich auf beliebiges $n$ -verallgemeinert werden. -Allerdings ist die Rechnung zur Bestimmung der Stützstellen und -Gewichte sehr mühsam. - -\subsubsection{Stützstellen und Orthogonalpolynome} -Sei $R_n=\{p(X)\in\mathbb{R}[X] \mid \deg p\le n\}$ der Vektorraum -der Polynome vom Grad $n$. - -\begin{satz} -\label{buch:integral:satz:gaussquadratur} -Sei $p$ ein Polynom vom Grad $n$, welches auf allen Polynomen in $R_{n-1}$ -orthogonal sind. -Seien ausserdem $x_0<x_1<\dots<x_n$ Stützstellen im Intervall $[-1,1]$ -und $A_i\in\mathbb{R}$ Gewichte derart dass -\[ -\int_{-1}^1 f(x)\,dx = -\sum_{i=0}^n A_if(x_i) -\] -für jedes Polynom $f$ vom Grad höchstens $2n-1$, dann sind die Zahlen -$x_i$ die Nullstellen des Polynoms $p$. -\end{satz} - -\begin{proof}[Beweis] -Sei $f(x)$ ein beliebiges Polynom vom Grad $2n-1$. -Nach dem Polynomdivisionsalgorithmus gibt es -Polynome $q,r\in R_{n-1}$ derart, dass $f=qp+r$. -Dann ist das Integral von $f$ gegeben durch -\[ -\int_{-1}^1 f(x)\,dx -= -\int_{-1}^1q(x) p(x)\,dx + \int_{-1}^1 r(x)\,dx -= -\langle q,p\rangle + \int_{-1}^1 r(x)\,dx. -\] -Da $p\perp R_{n-1}$ folgt insbesondere, dass $\langle q,p\rangle=0$. - -Da die Integrale auch aus den Werten in den Stützstellen berechnet -werden können, muss auch -\[ -0 -= -\int_{-1}^1 q(x)p(x)\,dx -= -\sum_{i=0}^n q(x_i)p(x_i) -\] -für jedes beliebige Polynom $q\in R_{n-1}$ gelten. -Da man für $q$ die Interpolationspolynome $l_j(x)$ verwenden -kann, den Grad $n-1$ haben, folgt -\[ -0 -= -\sum_{i=0}^n -l_j(x_i)p(x_i) -= -\sum_{i=0}^n \delta_{ij}p(x_i), -\] -die Stützstellen $x_i$ müssen also die Nullstellen des Polynoms -$p(x)$ sein. -\end{proof} - -Der Satz~\ref{buch:integral:satz:gaussquadratur} begründet das -{\em Gausssche Quadraturverfahren}. -Die in Abschnitt~\ref{buch:integral:subsection:orthogonale-polynome} -bestimmten Legendre-Polynome $P_n$ haben die im Satz -verlangte Eigenschaft, -dass sie auf allen Polynomen geringeren Grades orthogonal sind. -Wählt man die $n$ Nullstellen von $P_n$ als Stützstellen, erhält man -automatisch ein Integrationsverfahren, welches für Polynome vom Grad -$2n-1$ exakt ist. - -\begin{beispiel} -Das Legendre-Polynom $P_2(x) = \frac12(3x^2-1)$ hat die -Nullstellen $x=\pm1/\sqrt{3}$, dies sind genau die im Beispiel -auf Seite~\pageref{buch:integral:beispiel:gaussquadraturn1} befundenen -Sützstellen. -\end{beispiel} - -\subsubsection{Fehler der Gauss-Quadratur} -Das Gausssche Quadraturverfahren mit $n$ Stützstellen berechnet -Integrale von Polynomen bis zum Grad $2n-1$ exakt. -Für eine beliebige Funktion kann man die folgende Fehlerabschätzung -angeben \cite[theorem 7.3.4, p.~497]{numal}. - -\begin{satz} -Seien $x_i$ die Stützstellen und $A_i$ die Gewichte einer -Gaussschen Quadraturformel mit $n+1$ Stützstellen und sei $f$ -eine auf dem Interval $[-1,1]$ $2n+2$-mal stetig differenzierbare -Funktion, dann ist der $E$ Fehler des Integrals -\[ -\int_{-1}^1 f(x)\,dx = \sum_{i=0}^n A_i f(x_i) + E -\] -gegeben durch -\begin{equation} -E = \frac{f^{(2n+2)}}{(2n+2)!}\int_{-1}^1 l(x)^2\,dx, -\label{buch:integral:gaussquadratur:eqn:fehlerformel} -\end{equation} -wobei $l(x)=(x-x_0)(x-x_1)\dots(x-x_n)$ ist. -\end{satz} - -Dank dem Faktor $(2n+2)!$ im Nenner von -\eqref{buch:integral:gaussquadratur:eqn:fehlerformel} -geht der Fehler für grosses $n$ sehr schnell gegen $0$. -Man kann auch zeigen, dass die mit Gauss-Quadratur mit $n+1$ -Stützstellen berechneten Näherungswerte eines Integrals einer -stetigen Funktion $f(x)$ für $n\to\infty$ immer gegen den wahren -Wert des Integrals konvergieren. - -\subsubsection{Skalarprodukte mit Gewichtsfunktion} - +%Das Interpolationspolynom $p(x)$ lässt sich abera uch direkt +%angeben. +%Dazu konstruiert man zuerst die Polynome +%\[ +%l_i(x) +%= +%\frac{ +%(x-x_0)(x-x_1)\cdots\widehat{(x-x_i)}\cdots (x-x_n) +%}{ +%(x_i-x_0)(x_i-x_1)\cdots\widehat{(x_i-x_i)}\cdots (x_i-x_n) +%} +%\] +%vom Grad $n$, wobei der Hut bedeutet, dass diese Faktoren +%im Produkt wegzulassen sind. +%Die Polynome $l_i(x)$ haben die Eigenschaft +%\[ +%l_i(x_j) = \delta_{ij} +%= +%\begin{cases} +%1&\qquad i=j\\ +%0&\qquad\text{sonst}. +%\end{cases} +%\] +%Die Linearkombination +%\[ +%p(x) = \sum_{i=0}^n f(x_i)l_i(x) +%\] +%ist dann ein Polynom vom Grad $n$, welches am den Stellen $x_j$ +%die Werte +%\[ +%p(x_j) +%= +%\sum_{i=0}^n f(x_i)l_i(x_j) +%= +%\sum_{i=0}^n f(x_i)\delta_{ij} +%= +%f(x_j) +%\] +%hat, das Polynome $p(x)$ ist also das gesuchte Interpolationspolynom. +% +%\subsubsection{Integrationsverfahren auf der Basis von Interpolation} +%Das Integral einer stetigen Funktion $f(x)$ auf dem Intervall $[-1,1]$ +%kann mit Hilfe des Interpolationspolynoms approximiert werden. +%Wenn $|f(x)-p(x)|<\varepsilon$ ist im Intervall $[-1,1]$, dann gilt +%für die Integrale +%\[ +%\biggl|\int_{-1}^1 f(x)\,dx -\int_{-1}^1p(x)\,dx\biggr| +%\le +%\int_{-1}^1 |f(x)-p(x)|\,dx +%\le +%2\varepsilon. +%\] +%Ein Interpolationspolynom mit kleinem Fehler liefert also auch +%eine gute Approximation für das Integral. +% +%Da das Interpolationspolynome durch die Funktionswerte $f(x_i)$ +%bestimmt ist, muss auch das Integral allein aus diesen Funktionswerten +%berechnet werden können. +%Tatsächlich ist +%\begin{equation} +%\int_{-1}^1 p(x)\,dx +%= +%\int_{-1}^1 \sum_{i=0}^n f(x_i)l_i(x)\,dx +%= +%\sum_{i=0}^n f(x_i) +%\underbrace{\int_{-1}^1 +%l_i(x)\,dx}_{\displaystyle = A_i}. +%\label{buch:integral:gaussquadratur:eqn:Aidef} +%\end{equation} +%Das Integral von $f(x)$ wird also durch eine mit den Zahlen $A_i$ +%gewichtete Summe +%\[ +%\int_{-1}^1 f(x)\,dx +%\approx +%\sum_{i=1}^n f(x_i)A_i +%\] +%approximiert. +% +%\subsubsection{Integrationsverfahren, die für Polynome exakt sind} +%Ein Polynom vom Grad $2n$ hat $2n+1$ Koeffizienten. +%Um das Polynom durch ein Interpolationspolynom exakt wiederzugeben, +%braucht man $2n+1$ Stützstellen. +%Andererseits gilt +%\[ +%\int_{-1}^1 a_{2n}x^{2n} + a_{2n-1}x^{2n-1} + \dots + a_2x^2 + a_1x a_0\,dx +%= +%\int_{-1}^1 a_{2n}x^{2n} + a_{2n-2}x^{2n-2}+\dots +a_2x^2 +a_0\,dx, +%\] +%das Integral ist also bereits durch die $n+1$ Koeffizienten mit geradem +%Index bestimmt. +%Es sollte daher möglich sein, aus $n+1$ Funktionswerten eines beliebigen +%Polynoms vom Grad höchstens $2n$ an geeignet gewählten Stützstellen das +%Integral exakt zu bestimmen. +% +%\begin{beispiel} +%Wir versuchen dies für quadratische Polynome durchzuführen, also +%für $n=1$. +%Gesucht sind also zwei Werte $x_i$, $i=0,1$ und Gewichte $A_i$, $i=0,1$ +%derart, dass für jedes quadratische Polynome $p(x)=a_2x^2+a_1x+a_0$ +%das Integral durch +%\[ +%\int_{-1}^1 p(x)\,dx +%= +%A_0 p(x_0) + A_1 p(x_1) +%\] +%gebeben ist. +%Indem wir für $p(x)$ die Polynome $1$, $x$, $x^2$ und $x^3$ einsetzen, +%erhalten wir vier Gleichungen +%\[ +%\begin{aligned} +%p(x)&=\rlap{$1$}\phantom{x^2}\colon& 2 &= A_0\phantom{x_0}+ A_1 \\ +%p(x)&=x^{\phantom{2}}\colon& 0 &= A_0x_0 + A_1x_1 \\ +%p(x)&=x^2\colon& \frac23 &= A_0x_0^2 + A_1x_1^2\\ +%p(x)&=x^3\colon& 0 &= A_0x_0^3 + A_1x_1^3. +%\end{aligned} +%\] +%Dividiert man die zweite und vierte Gleichung in der Form +%\[ +%\left. +%\begin{aligned} +%A_0x_0 &= -A_1x_1\\ +%A_0x_0^2 &= -A_1x_1^2 +%\end{aligned} +%\quad +%\right\} +%\quad +%\Rightarrow +%\quad +%x_0^2=x_1^2 +%\quad +%\Rightarrow +%\quad +%x_1=-x_0. +%\] +%Indem wir dies in die zweite Gleichung einsetzen, finden wir +%\[ +%0 = A_0x_0 + A_1x_1 = A_0x_1 -A_1x_0 = (A_0-A_1)x_0 +%\quad\Rightarrow\quad +%A_0=A_1. +%\] +%Aus der ersten Gleichung folgt jetzt +%\[ +%2= A_0+A_1 = 2A_0 \quad\Rightarrow\quad A_0 = 1. +%\] +%Damit bleiben nur noch die Werte von $x_i$ zu bestimmen, was +%mit Hilfe der zweiten Gleichung geschehen kann: +%\[ +%\frac23 = A_0x_0^2 + A_1x_1^2 = 2x_0^2 +%\quad\Rightarrow\quad +%x_0 = \frac{1}{\sqrt{3}}, x_1 = -\frac{1}{\sqrt{3}} +%\] +%Damit ist das Problem gelöst: das Integral eines Polynoms vom Grad 3 +%im Interval $[-1,1]$ ist exakt gegeben durch +%\[ +%\int_{-1}^1 p(x)\,dx +%= +%p\biggl(-\frac{1}{\sqrt{3}}\biggr) +%+ +%p\biggl(\frac{1}{\sqrt{3}}\biggr). +%\] +%Das Integral kann also durch nur zwei Auswertungen des Polynoms +%exakt bestimmt werden. +% +%Im Laufe der Lösung des Gleichungssystems wurden die Gewichte $A_i$ +%mit bestimmt. +%Es ist aber auch möglich, die Gewichte zu bestimmen, wenn man die +%Stützstellen kennt. +%Nach \eqref{buch:integral:gaussquadratur:eqn:Aidef} +%sind sie die $A_i$ gegeben als Integrale der Polynome +%$l_i(x)$, die im vorliegenden Fall linear sind: +%\begin{align*} +%l_0(x) +%&= +%\frac{x-x_1}{x_0-x_1} +%= +%\frac{x-\frac1{\sqrt{3}}}{-\frac{2}{\sqrt{3}}} +%= +%\frac12(1-\sqrt{3}x) +%\\ +%l_1(x) +%&= +%\frac{x-x_0}{x_1-x_0} +%= +%\frac{x+\frac1{\sqrt{3}}}{\frac{2}{\sqrt{3}}} +%= +%\frac12(1+\sqrt{3}x) +%\end{align*} +%Diese haben die Integrale +%\[ +%\int_{-1}^1\frac12(1\pm\sqrt{3}x)\,dx +%= +%\int_{-1}^1 \frac12\,dx +%= +%1, +%\] +%da das Polynom $x$ verschwindendes Integral hat. +%Dies stimmt mit $A_0=A_1=1$ überein. +%\label{buch:integral:beispiel:gaussquadraturn1} +%\end{beispiel} +% +%Das eben vorgestellt Verfahren kann natürlich auf beliebiges $n$ +%verallgemeinert werden. +%Allerdings ist die Rechnung zur Bestimmung der Stützstellen und +%Gewichte sehr mühsam. +% +%\subsubsection{Stützstellen und Orthogonalpolynome} +%Sei $R_n=\{p(X)\in\mathbb{R}[X] \mid \deg p\le n\}$ der Vektorraum +%der Polynome vom Grad $n$. +% +%\begin{satz} +%\label{buch:integral:satz:gaussquadratur} +%Sei $p$ ein Polynom vom Grad $n$, welches auf allen Polynomen in $R_{n-1}$ +%orthogonal sind. +%Seien ausserdem $x_0<x_1<\dots<x_n$ Stützstellen im Intervall $[-1,1]$ +%und $A_i\in\mathbb{R}$ Gewichte derart dass +%\[ +%\int_{-1}^1 f(x)\,dx = +%\sum_{i=0}^n A_if(x_i) +%\] +%für jedes Polynom $f$ vom Grad höchstens $2n-1$, dann sind die Zahlen +%$x_i$ die Nullstellen des Polynoms $p$. +%\end{satz} +% +%\begin{proof}[Beweis] +%Sei $f(x)$ ein beliebiges Polynom vom Grad $2n-1$. +%Nach dem Polynomdivisionsalgorithmus gibt es +%Polynome $q,r\in R_{n-1}$ derart, dass $f=qp+r$. +%Dann ist das Integral von $f$ gegeben durch +%\[ +%\int_{-1}^1 f(x)\,dx +%= +%\int_{-1}^1q(x) p(x)\,dx + \int_{-1}^1 r(x)\,dx +%= +%\langle q,p\rangle + \int_{-1}^1 r(x)\,dx. +%\] +%Da $p\perp R_{n-1}$ folgt insbesondere, dass $\langle q,p\rangle=0$. +% +%Da die Integrale auch aus den Werten in den Stützstellen berechnet +%werden können, muss auch +%\[ +%0 +%= +%\int_{-1}^1 q(x)p(x)\,dx +%= +%\sum_{i=0}^n q(x_i)p(x_i) +%\] +%für jedes beliebige Polynom $q\in R_{n-1}$ gelten. +%Da man für $q$ die Interpolationspolynome $l_j(x)$ verwenden +%kann, den Grad $n-1$ haben, folgt +%\[ +%0 +%= +%\sum_{i=0}^n +%l_j(x_i)p(x_i) +%= +%\sum_{i=0}^n \delta_{ij}p(x_i), +%\] +%die Stützstellen $x_i$ müssen also die Nullstellen des Polynoms +%$p(x)$ sein. +%\end{proof} +% +%Der Satz~\ref{buch:integral:satz:gaussquadratur} begründet das +%{\em Gausssche Quadraturverfahren}. +%Die in Abschnitt~\ref{buch:integral:subsection:orthogonale-polynome} +%bestimmten Legendre-Polynome $P_n$ haben die im Satz +%verlangte Eigenschaft, +%dass sie auf allen Polynomen geringeren Grades orthogonal sind. +%Wählt man die $n$ Nullstellen von $P_n$ als Stützstellen, erhält man +%automatisch ein Integrationsverfahren, welches für Polynome vom Grad +%$2n-1$ exakt ist. +% +%\begin{beispiel} +%Das Legendre-Polynom $P_2(x) = \frac12(3x^2-1)$ hat die +%Nullstellen $x=\pm1/\sqrt{3}$, dies sind genau die im Beispiel +%auf Seite~\pageref{buch:integral:beispiel:gaussquadraturn1} befundenen +%Sützstellen. +%\end{beispiel} +% +%\subsubsection{Fehler der Gauss-Quadratur} +%Das Gausssche Quadraturverfahren mit $n$ Stützstellen berechnet +%Integrale von Polynomen bis zum Grad $2n-1$ exakt. +%Für eine beliebige Funktion kann man die folgende Fehlerabschätzung +%angeben \cite[theorem 7.3.4, p.~497]{buch:numal}. +% +%\begin{satz} +%Seien $x_i$ die Stützstellen und $A_i$ die Gewichte einer +%Gaussschen Quadraturformel mit $n+1$ Stützstellen und sei $f$ +%eine auf dem Interval $[-1,1]$ $2n+2$-mal stetig differenzierbare +%Funktion, dann ist der $E$ Fehler des Integrals +%\[ +%\int_{-1}^1 f(x)\,dx = \sum_{i=0}^n A_i f(x_i) + E +%\] +%gegeben durch +%\begin{equation} +%E = \frac{f^{(2n+2)}(\xi)}{(2n+2)!}\int_{-1}^1 l(x)^2\,dx, +%\label{buch:integral:gaussquadratur:eqn:fehlerformel} +%\end{equation} +%wobei $l(x)=(x-x_0)(x-x_1)\dots(x-x_n)$ und $\xi$ ein geeigneter +%Wert im Intervall $[-1,1]$ ist. +%\end{satz} +% +%Dank dem Faktor $(2n+2)!$ im Nenner von +%\eqref{buch:integral:gaussquadratur:eqn:fehlerformel} +%geht der Fehler für grosses $n$ sehr schnell gegen $0$. +%Man kann auch zeigen, dass die mit Gauss-Quadratur mit $n+1$ +%Stützstellen berechneten Näherungswerte eines Integrals einer +%stetigen Funktion $f(x)$ für $n\to\infty$ immer gegen den wahren +%Wert des Integrals konvergieren. +% +%\begin{table} +%\def\u#1{\underline{#1}} +%\centering +%\begin{tabular}{|>{$}c<{$}|>{$}r<{$}|>{$}r<{$}|} +%\hline +% n & \text{Gauss-Quadratur} & \text{Trapezregel} \\ +%\hline +%\phantom{0}2 & 0.\u{95}74271077563381 & 0.\u{95}63709682242596 \\ +%\phantom{0}4 & 0.\u{95661}28333449730 & 0.\u{956}5513401768598 \\ +%\phantom{0}6 & 0.\u{9566114}812034364 & 0.\u{956}5847489712136 \\ +%\phantom{0}8 & 0.\u{956611477}5028123 & 0.\u{956}5964425360520 \\ +% 10 & 0.\u{9566114774905}637 & 0.\u{9566}018550715587 \\ +% 12 & 0.\u{956611477490518}7 & 0.\u{9566}047952369826 \\ +% 14 & 0.\u{95661147749051}72 & 0.\u{9566}065680717177 \\ +% 16 & 0.\u{956611477490518}7 & 0.\u{9566}077187127541 \\ +% 18 & 0.\u{956611477490518}3 & 0.\u{9566}085075898731 \\ +% 20 & 0.\u{956611477490518}4 & 0.\u{9566}090718697414 \\ +%\hline +% \infty & 0.9566114774905183 & 0.9566114774905183 \\ +%\hline +%\end{tabular} +%\caption{Integral von $\sqrt{1-x^2}$ zwischen $-\frac12$ und $\frac12$ +%berechnet mit Gauss-Quadratur und der Trapezregel, aber mit zehnmal +%so vielen Stützstellen. +%Bereits mit 12 Stützstellen erreicht die Gauss-Quadratur +%Maschinengenauigkeit, die Trapezregel liefert auch mit 200 Stützstellen +%nicht mehr als 4 korrekte Nachkommastellen. +%\label{buch:integral:gaussquadratur:table0.5}} +%\end{table} +% +%%\begin{table} +%%\def\u#1{\underline{#1}} +%%\centering +%%\begin{tabular}{|>{$}c<{$}|>{$}r<{$}|>{$}r<{$}|} +%%\hline +%% n & \text{Gauss-Quadratur} & \text{Trapezregel} \\ +%%\hline +%%\phantom{0}2 & 1.\u{5}379206741571556 & 1.\u{5}093105464758343 \\ +%%\phantom{0}4 & 1.\u{51}32373472933831 & 1.\u{51}13754509594814 \\ +%%\phantom{0}6 & 1.\u{512}1624557410367 & 1.\u{51}17610879524799 \\ +%%\phantom{0}8 & 1.\u{51207}93479994321 & 1.\u{51}18963282632112 \\ +%% 10 & 1.\u{51207}13859966004 & 1.\u{51}19589735776959 \\ +%% 12 & 1.\u{512070}5317779943 & 1.\u{51}19930161260693 \\ +%% 14 & 1.\u{5120704}334802813 & 1.\u{5120}135471596636 \\ +%% 16 & 1.\u{5120704}216176006 & 1.\u{5120}268743889558 \\ +%% 18 & 1.\u{5120704}201359081 & 1.\u{5120}360123137213 \\ +%% 20 & 1.\u{5120704199}459651 & 1.\u{5120}425490275837 \\ +%%\hline +%% \infty & 1.5120704199172947 & 1.5120704199172947 \\ +%%\hline +%%\end{tabular} +%%\end{table} +% +%%\begin{table} +%%\def\u#1{\underline{#1}} +%%\centering +%%\begin{tabular}{|>{$}c<{$}|>{$}r<{$}|>{$}r<{$}|} +%%\hline +%% n & \text{Gauss-Quadratur} & \text{Trapezregel} \\ +%%\hline +%%\phantom{0}2 & 1.\u{}6246862220133462 & 1.\u{5}597986803933712 \\ +%%\phantom{0}4 & 1.\u{5}759105515463101 & 1.\u{56}63563456168101 \\ +%%\phantom{0}6 & 1.\u{5}706630058381434 & 1.\u{56}77252866190838 \\ +%%\phantom{0}8 & 1.\u{56}94851106536780 & 1.\u{568}2298707696152 \\ +%% 10 & 1.\u{56}91283195332679 & 1.\u{568}4701957758742 \\ +%% 12 & 1.\u{56}90013806299465 & 1.\u{568}6030805941198 \\ +%% 14 & 1.\u{5689}515434853885 & 1.\u{568}6841603070025 \\ +%% 16 & 1.\u{5689}306507843050 & 1.\u{568}7372230731711 \\ +%% 18 & 1.\u{5689}214761291217 & 1.\u{568}7738235496322 \\ +%% 20 & 1.\u{56891}73062385982 & 1.\u{568}8001228530786 \\ +%%\hline +%% \infty & 1.5689135396691616 & 1.5689135396691616 \\ +%%\hline +%%\end{tabular} +%%\end{table} +% +%\begin{table} +%\def\u#1{\underline{#1}} +%\centering +%\begin{tabular}{|>{$}c<{$}|>{$}r<{$}|>{$}r<{$}|} +%\hline +% n & \text{Gauss-Quadratur} & \text{Trapezregel} \\ +%\hline +%\phantom{0}2 & 1.\u{}6321752373234928 & 1.\u{5}561048774629949 \\ +%\phantom{0}4 & 1.\u{57}98691557134743 & 1.\u{5}660124134617943 \\ +%\phantom{0}6 & 1.\u{57}35853681692993 & 1.\u{5}683353001877542 \\ +%\phantom{0}8 & 1.\u{57}19413565928206 & 1.\u{5}692627503425400 \\ +% 10 & 1.\u{57}13388119633434 & 1.\u{5}697323578543481 \\ +% 12 & 1.\u{57}10710489948883 & 1.\u{570}0051217458713 \\ +% 14 & 1.\u{570}9362135398341 & 1.\u{570}1784766276063 \\ +% 16 & 1.\u{570}8621102742815 & 1.\u{570}2959121005231 \\ +% 18 & 1.\u{570}8186779483588 & 1.\u{570}3793521168343 \\ +% 20 & 1.\u{5707}919411931615 & 1.\u{570}4408749735932 \\ +%\hline +% \infty & 1.5707367072605671 & 1.5707367072605671 \\ +%\hline +%\end{tabular} +%\caption{Integral von $\sqrt{1-x^2}$ zwischen $-0.999$ und $0.999$ +%berechnet mit Gauss-Quadratur und der Trapezregel, aber mit zehnmal +%so vielen Stützstellen. +%Wegen der divergierenden Steigung des Integranden bei $\pm 1$ tun +%sich beide Verfahren sehr schwer. +%Trotzdem erreich die Gauss-Quadrator 4 korrekte Nachkommastellen +%mit 20 Stütztstellen, während die Trapezregel auch mit 200 Stützstellen +%nur 3 korrekte Nachkommastellen findet. +%\label{buch:integral:gaussquadratur:table0.999}} +%\end{table} +% +%\begin{figure} +%\centering +%\includegraphics{chapters/060-integral/gq/gq.pdf} +%\caption{Approximationsfehler des +%Integrals~\eqref{buch:integral:gaussquadratur:bspintegral} +%in Abhängigkeit von $a$. +%Die Divergenz der Ableitung des Integranden an den Intervallenden +%$\pm 1$ führt zu schlechter Konvergenz des Verfahrens, wenn $a$ +%nahe an $1$ ist. +%\label{buch:integral:gaussquadratur:fehler}} +%\end{figure} +% +%Zur Illustration der Genauigkeit der Gauss-Quadratur berechnen wir +%das Integral +%\begin{equation} +%\int_{-a}^a \sqrt{1-x^2}\,dx +%= +%\arcsin a + a \sqrt{1-a^2} +%\label{buch:integral:gaussquadratur:bspintegral} +%\end{equation} +%mit Gauss-Quadratur einerseits und dem Trapezverfahren +%andererseits. +%Da Gauss-Quadratur mit sehr viel weniger Sützstellen auskommt, +%berechnen wir die Trapeznäherung mit zehnmal so vielen Stützstelln. +%In den Tabellen~\ref{buch:integral:gaussquadratur:table0.5} +%und +%\ref{buch:integral:gaussquadratur:table0.999} +%sind die Resultate zusammengestellt. +%Für $a =\frac12$ zeigt +%Tabelle~\ref{buch:integral:gaussquadratur:table0.5} +%die sehr schnelle Konvergenz der Gauss-Quadratur, schon mit +%12 Stützstellen wird Maschinengenauigkeit erreicht. +%Das Trapezverfahren dagegen erreicht auch mit 200 Stützstellen nur +%4 korrekte Nachkommastellen. +% +%An den Stellen $x=\pm 1$ divergiert die Ableitung des Integranden +%des Integrals \eqref{buch:integral:gaussquadratur:bspintegral}. +%Da grösste und kleinste Stützstelle der Gauss-Quadratur immer +%deutlich vom Rand des Intervalls entfernt ist, kann das Verfahren +%diese ``schwierigen'' Stellen nicht erkennen. +%Tabelle~\ref{buch:integral:gaussquadratur:table0.999} zeigt, wie +%die Konvergenz des Verfahrens in diesem Fall sehr viel schlechter ist. +%Dies zeigt auch der Graph in +%Abbildung~\ref{buch:integral:gaussquadratur:fehler}. +% +%\subsubsection{Skalarprodukte mit Gewichtsfunktion} +\input{chapters/060-integral/gaussquadratur.tex} diff --git a/buch/chapters/references.bib b/buch/chapters/references.bib index 20215a9..80a8379 100644 --- a/buch/chapters/references.bib +++ b/buch/chapters/references.bib @@ -78,7 +78,7 @@ @book{buch:ab, title = {$A = B$}, - author = {Marko Petkov\v sek, Herbert S. Wilf and Doron Zeilberger}, + author = {Marko Petkov\v{s}ek and Herbert S. Wilf and Doron Zeilberger}, year = 1996, publisher = {A K Peters, Ltd}, ISBN = {978-1-56881-063-8} |