From 0c40af1b99a0b0f60be8786b65c277ce7813ee12 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20M=C3=BCller?= Date: Sun, 28 Nov 2021 19:52:52 +0100 Subject: gauss quadratur --- buch/chapters/060-integral/fehlerfunktion.tex | 52 ++ buch/chapters/060-integral/orthogonal.tex | 826 +++++++++++++++++++++++++- buch/chapters/references.bib | 11 + 3 files changed, 887 insertions(+), 2 deletions(-) diff --git a/buch/chapters/060-integral/fehlerfunktion.tex b/buch/chapters/060-integral/fehlerfunktion.tex index 37a007d..397615e 100644 --- a/buch/chapters/060-integral/fehlerfunktion.tex +++ b/buch/chapters/060-integral/fehlerfunktion.tex @@ -576,6 +576,58 @@ Diese Reihenentwicklung ist sehr effizient für kleine Werte von $x$. Für grosse Werte von $x$ entstehen aber sehr grosse Zwischenterme in der Reihe, was zu Auslöschung und damit zu Genauigkeitsverlust führt. +\subsubsection{Hypergeometrische Funktion} +Die Taylor-Reihe~\eqref{buch:integrale:eqn:erfreihe} der Fehlerfunktion +kann auch mit Hilfe hypergeometrischer Funktionen geschrieben werden. +Da nur ungerade Potenzen vorkommen, klammern wir zunächst einen gemeinsamen +Faktor $x$ aus: +\[ +\operatorname{erf}(x) += +\frac{2x}{\sqrt{\pi}} +\sum_{k=0}^\infty +\frac{1}{2k+1} +\frac{(-x^2)^k}{k!}. +\] +Der Bruch $1/(2k+1)$ muss jetzt noch mit Hilfe von Pochhammer-Symbolen +geschrieben werden. +Dazu beachten wir, dass +\begin{align*} +\frac{1}{2k+1} +&= +\frac12 +\frac{1}{\frac32+k-1} +\\ +&= +\frac12 +\frac{ +\frac32(\frac32+1)(\frac32+2)\dots(\frac32+k-2)\phantom{(\frac32+k-1)} +}{ +\frac32(\frac32+1)(\frac32+2)\dots(\frac32+k-2)(\frac32+k-1) +} +\\ +&= +\frac{ +\frac12(\frac12+1)(\frac12+2)(\frac12+3)\dots(\frac12+k-1) +}{ +\frac32(\frac32+1)(\frac32+2)\dots(\frac32+k-2)(\frac32+k-1) +} +\\ +&= +\frac{(\frac12)_k}{(\frac32)_k}. +\end{align*} +Somit ist die Fehlerfunktion als hypergeometrische Funktion +\[ +\operatorname{erf}(x) += +\frac{2x}{\sqrt{\pi}}\sum_{k=0}^\infty +\frac{(\frac12)_k}{(\frac32)_k}\frac{(-x^2)^k}{k!} += +\frac{2x}{\sqrt{\pi}}\, +\mathstrut_1F_1({\textstyle\frac12};{\textstyle\frac32};-x^2). +\] +gegeben. + \subsubsection{Kettenbruchentwicklung} Besonders für grosse $x$ interessiert man sich mehr für $\operatorname{erfc}(x)$ als für $\operatorname{erf}(x)$. diff --git a/buch/chapters/060-integral/orthogonal.tex b/buch/chapters/060-integral/orthogonal.tex index 6300e26..4ffbbde 100644 --- a/buch/chapters/060-integral/orthogonal.tex +++ b/buch/chapters/060-integral/orthogonal.tex @@ -3,7 +3,7 @@ % % (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule % -\section{Orthogonale Polynome +\section{Orthogonalität \label{buch:integral:section:orthogonale-polynome}} Die Fourier-Theorie basiert auf der Idee, Funktionen durch Funktionenreihen mit Summanden zu bilden, die im Sinne eines @@ -14,12 +14,834 @@ Differentialgleichungen. Besonders interessant wird die Situation, wenn die Funktionen Polynome sind. +% +% Skalarprodukt +% \subsection{Skalarprodukt} -\subsection{Definition} +Der reelle Vektorraum $\mathbb{R}^n$ trägt das Skalarprodukt +\[ +\langle\;,\;\rangle +\colon +\mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R} +: +(x,y)\mapsto \langle x, y\rangle = \sum_{k=1}^n x_iy_k, +\] +welches viele interessante Anwendungen ermöglicht. +Eine orthonormierte Basis macht es zum Beispiel besonders leicht, +eine Zerlegung eines Vektors in dieser Basis zu finden. +In diesem Abschnitt soll zunächst an die Eigenschaften erinnert +werden, die zu einem nützlichen + +\subsubsection{Eigenschaften eines Skalarproduktes} +Das Skalarprodukt erlaubt auch, die Länge eines Vektors $v$ +als $|v| = \sqrt{\langle v,v\rangle}$ zu definieren. +Dies funktioniert natürlich nur, wenn die Wurzel auch immer +definiert ist, d.~h.~das Skalarprodukt eines Vektors mit sich +selbst darf nicht negativ sein. +Dazu dient die folgende Definition. + +\begin{definition} +Sei $V$ ein reeller Vektorraum. +Eine bilineare Abbildung +\[ +\langle\;,\;\rangle +\colon +V\times V +\to +\mathbb{R} +: +(u,v) \mapsto \langle u,v\rangle. +\] +heisst {\em positiv definit}, wenn für alle Vektoren $v \in V$ mit +$v\ne 0 \Rightarrow \langle v,v\rangle > 0$ +Die {\em Norm} eines Vektors $v$ ist +$|v|=\sqrt{\langle v,v\rangle}$. +\end{definition} + +Damit man mit dem Skalarprodukt sinnvoll rechnen kann, ist ausserdem +erforderlich, dass es eine einfache Beziehung zwischen +$\langle x,y\rangle$ und $\langle y,x\rangle$ gibt. + +\begin{definition} +Ein {\em Skalarprodukt} auf einem reellen Vektorraum $V$ ist eine +positiv definite, symmetrische bilineare Abbildung +\[ +\langle\;,\;\rangle +\colon +V\times V +\to +\mathbb{R} +: +(u,v) \mapsto \langle u,v\rangle. +\] +\end{definition} + +Das Skalarprodukt $\langle u,v\rangle=u^tv$ auf dem Vektorraum +$\mathbb{R}^n$ erfüllt die Definition ganz offensichtlich, +sie führt auf die Komponentendarstellung +\[ +\langle u,v\rangle = u^tv = \sum_{k=1}^n u_iv_i. +\] +Weitere Skalarprodukte ergeben ergeben sich mit jeder symmetrischen, +positiv definiten Matrix $G$ und der Definition +$\langle u,v\rangle_G=u^tGv$. +Ein einfacher Spezialfall tritt auf, wenn $G$ eine Diagonalmatrix +$\operatorname{diag}(w_1,\dots,w_n)$ +mit positiven Einträgen $w_i>0$ auf der Diagonalen ist. +In diesem Fall schreiben wir +\[ +\langle u,v\rangle_w += +u^t\operatorname{diag}(w_1,\dots,w_n)v += +\sum_{k=1}^n u_iv_i\,w_i +\] +und nennen $\langle \;,\;\rangle_w$ das {\em gewichtete Skalarprodukt} +mit {\em Gewichten $w_i$}. + +\subsubsection{Skalarprodukte auf Funktionenräumen} +Das Integral ermöglicht jetzt, ein Skalarprodukt auf dem reellen +Vektorraum der stetigen Funktionen auf einem Intervall zu definieren. + +\begin{definition} +Sei $V$ der reelle Vektorraum $C([a,b])$ der reellwertigen, stetigen +Funktion auf dem Intervall $[a,b]$. +Dann ist +\[ +\langle\;,\;\rangle +\colon +C([a,b]) \times C([a,b]) \to \mathbb{R} +: +(f,g) \mapsto \langle f,g\rangle = \int_a^b f(x)g(x)\,dx. +\] +ein Skalarprodukt. +\end{definition} + +Die Definition ist offensichtlich symmetrisch in $f$ und $g$ und +aus den Eigenschaften des Integrals ist klar, dass das Produkt +bilinear ist: +\begin{align*} +\langle \lambda_1 f_1+\lambda_2f_2,g\rangle +&= +\int_a^b (\lambda_1f_(x) +\lambda_2f_2(x))g(x)\,dx += +\lambda_1\int_a^b f_1(x) g(x)\,dx ++ +\lambda_2\int_a^b f_2(x) g(x)\,dx +\\ +&= +\lambda_1\langle f_1,g\rangle ++ +\lambda_2\langle f_2,g\rangle. +\end{align*} +Ausserdem ist es positiv definit, denn wenn $f(x_0) \ne 0$ ist, +dann gibt es wegen der Stetigkeit von $f$ eine Umgebung +$U=[x_0-\varepsilon,x_0+\varepsilon]$, derart, dass $|f(x)| > \frac12|f(x_0)|$ +ist für alle $x\in U$. +Somit ist das Integral +\[ +\langle f,f\rangle += +\int_a^b |f(x)|^2\,dx +\ge +\int_{x_0-\varepsilon}^{x_0+\varepsilon} |f(x)|^2\,dx +\ge +\int_{x_0-\varepsilon}^{x_0+\varepsilon} \frac14|f(x_0)|^2\,dx += +\frac{1}{4}|f(x_0)|^2\cdot 2\varepsilon += +\frac{|f(x_0)|\varepsilon}{2} +>0, +\] +was beweist, dass $\langle\;,\;\rangle$ positiv definit und damit +ein Skalarprodukt ist. + +Die Definition kann noch etwas verallgemeinert werden, indem +die Funktionswerte nicht überall auf dem Definitionsbereich +gleich gewichtet werden. + +\begin{definition} +Sei $w\colon [a,b]\to \mathbb{R}^+$ eine positive, stetige Funktion, +dann ist +\[ +\langle\;,\;\rangle_w +\colon +C([a,b]) \times C([a,b]) \to \mathbb{R} +: +(f,g) \mapsto \langle f,g\rangle_w = \int_a^b f(x)g(x)\,w(x)\,dx. +\] +das {\em gewichtete Skalarprodukt} mit {\em Gewichtsfunktion $w(x)$}. +\end{definition} + +\subsubsection{Gram-Schmidt-Orthonormalisierung} +In einem reellen Vektorraum $V$ mit Skalarprodukt $\langle\;\,,\;\rangle$ +kann aus einer beleibigen Basis $b_1,\dots,b_n$ mit Hilfe des +Gram-Schmidtschen Orthogonalisierungsverfahrens immer eine +orthonormierte Basis $\tilde{b}_1,\dots,\tilde{b}_n$ Basis +gewonnen werden. +Es stellt sicher, dass für alle $k\le n$ gilt +\[ +\langle b_1,\dots,b_k\rangle += +\langle \tilde{b}_1,\dots,\tilde{b}_k\rangle. +\] +Zur Vereinfachung der Formeln schreiben wir $v^0=v/|v|$ für einen zu +$v$ parallelen Einheitsvektor. +Die Vektoren $\tilde{b}_i$ können mit Hilfe der Formeln +\begin{align*} +\tilde{b}_1 +&= +(b_1)^0 +\\ +\tilde{b}_2 +&= +\bigl( +b_2 +- +\langle \tilde{b}_1,b_2\rangle \tilde{b}_1 +\bigr)^0 +\\ +\tilde{b}_3 +&= +\bigl( +b_3 +- +\langle \tilde{b}_1,b_3\rangle \tilde{b}_1 +- +\langle \tilde{b}_2,b_3\rangle \tilde{b}_2 +\bigr)^0 +\\ +&\;\vdots +\\ +\tilde{b}_n +&= +\bigl( +b_n +- +\langle \tilde{b}_1,b_n\rangle \tilde{b}_1 +- +\langle \tilde{b}_2,b_n\rangle \tilde{b}_2 +-\dots +- +\langle \tilde{b}_{n-1},b_n\rangle \tilde{b}_{n-1} +\bigr)^0 +\end{align*} +iterativ berechnet werden. +Dieses Verfahren lässt sich auch auf Funktionenräume anwenden. + +Die Normierung ist nicht unbedingt nötig und manchmal unangenehm, +da die Norm unschöne Quadratwurzeln einführt. +Falls es genügt, eine orthogonale Basis zu finden, kann darauf +verzichtet werden, bei der Orthogonalisierung muss aber berücksichtigt +werden, dass die Vektoren $\tilde{b}_i$ jetzt nicht mehr Einheitslänge +haben. +Die Formeln +\begin{align*} +\tilde{b}_0 +&= +b_0 +\\ +\tilde{b}_1 +&= +b_1 +- +\frac{\langle b_1,\tilde{b}_0\rangle}{\langle \tilde{b}_0,\tilde{b}_0\rangle}\tilde{b}_0 +\\ +\tilde{b}_2 +&= +b_2 +- +\frac{\langle b_2,\tilde{b}_0\rangle}{\langle \tilde{b}_0,\tilde{b}_0\rangle}\tilde{b}_0 +- +\frac{\langle b_2,\tilde{b}_1\rangle}{\langle \tilde{b}_1,\tilde{b}_1\rangle}\tilde{b}_1 +\\ +&\;\vdots +\\ +\tilde{b}_n +&= +b_n +- +\frac{\langle b_n,\tilde{b}_0\rangle}{\langle \tilde{b}_0,\tilde{b}_0\rangle}\tilde{b}_0 +- +\frac{\langle b_n,\tilde{b}_1\rangle}{\langle \tilde{b}_1,\tilde{b}_1\rangle}\tilde{b}_1 +- +\dots +- +\frac{\langle b_n,\tilde{b}_{n-1}\rangle}{\langle \tilde{b}_{n-1},\tilde{b}_{n-1}\rangle}\tilde{b}_{n-1}. +\end{align*} +berücksichtigen dies. + + +% +% Orthogonale Polynome +% +\subsection{Orthogonale Polynome +\label{buch:integral:subsection:orthogonale-polynome}} +Die Polynome $1,x,x^2,\dots,x^n$ bilden eine Basis des Vektorraums +der Polynome vom Grad $\le n$. +Bezüglich des Skalarproduktes +\[ +\langle p,q\rangle += +\int_{-1}^1 p(x)q(x)\,dx +\] +sind sie jedoch nicht orthogonal, denn es ist +\[ +\langle x^i,x^j\rangle += +\int_{-1}^1 x^{i+j}\,dx += +\biggl[\frac{x^{i+j+1}}{i+j+1}\biggr]_{-1}^1 += +\begin{cases} +\frac{2}{i+j+1}&\qquad\text{$i+j$ gerade}\\ + 0&\qquad\text{$i+j$ ungerade}. +\end{cases} +\] +Wir können daher das Gram-Schmidtsche Orthonormalisierungsverfahren +anwenden, um eine orthogonale Basis von Polynomen zu finden, was +wir im Folgenden tun wollen. + +% XXX Orthogonalisierungsproblem so formulieren, dass klar wird, +% XXX dass man ein "Normierungskriterium braucht. + +Da wir auf die Normierung verzichten, brauchen wir ein anderes +Kriterium, welches die Polynome eindeutig festlegen kann. +Wir bezeichnen das Polynom vom Grad $n$, das bei diesem Prozess +entsteht, mit $P_n(x)$ und legen willkürlich aber traditionskonform +fest, dass $P_n(1)=1$ sein soll. + +Das Skalarprodukt berechnet ein Integral eines Produktes von zwei +Polynomen über das symmetrische Interval $[-1,1]$. +Ist die eine gerade und die andere ungerade, dann ist das +Produkt eine ungerade Funktion und das Skalarprodukt verschwindet. +Sind beide Funktionen gerade oder ungerade, dann ist das Produkt +gerade und das Skalarprodukt ist im Allgmeinen von $0$ verschieden. +Dies zeigt, dass es tatsächlich etwas zu Orthogonalisieren gibt. + +Die ersten beiden Funktionen sind das konstante Polynom $1$ und +das Polynome $x$. +Nach obiger Beobachtung ist das Skalarprodukt $\langle 1,x\rangle=0$, +also ist $P_1(x)=x$. + +\begin{lemma} +Die Polynome $P_{2n}(x)$ sind gerade, die Polynome $P_{2n+1}(x)$ sind +ungerade Funktionen von $x$. +\end{lemma} + +\begin{proof}[Beweis] +Wir verwenden vollständige Induktion nach $n$. +Wir wissen bereits, dass $P_0(x)=1$ und $P_1(x)=x$ die verlangten +Symmetrieeigenschaften haben. +Im Sinne der Induktionsannahme nehmen wir daher an, dass die +Symmetrieeigenschaften für $P_k(x)$, $k{$}c<{$}|>{$}l<{$}|} +\hline +n&P_n(x)\\ +\hline + 0&1 +\\ + 1&x +\\ + 2&\frac12(3x^2-1) +\\ + 3&\frac12(5x^3-3x) +\\ + 4&\frac18(35x^4-30x^2+3) +\\ + 5&\frac18(63x^5-70x^3+15x) +\\ + 6&\frac1{16}(231x^6-315x^4+105x^2-5) +\\ + 7&\frac1{16}(429x^7-693x^5+315x^3-35x) +\\ + 8&\frac1{128}(6435x^8-12012x^6+6930x^4-1260x^2+35) +\\ + 9&\frac1{128}(12155x^9-25740x^7+18018x^5-4620x^3+315x) +\\ +10&\frac1{256}(46189x^{10}-109395x^8+90090x^6-30030x^4+3465x^2-63) +\\ +\hline +\end{tabular} +\caption{Die Legendre-Polynome $P_n(x)$ für $n=0,1,\dots,10$ sind +orthogonale Polynome vom Grad $n$, die den Wert $P_n(1)=1$ haben. +\label{buch:integral:table:legendre-polynome}} +\end{table} + +Die so konstruierten Polynome heissen die {\em Legendre-Polynome}. +Durch weitere Durchführung des Verfahrens liefert die Polynome in +Tabelle~\ref{buch:integral:table:legendre-polynome}. + + +% +% Differentialgleichungen +% \subsection{Orthogonale Polynome und Differentialgleichungen} \subsubsection{Legendre-Differentialgleichung} \subsubsection{Legendre-Polyome} \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. + +\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