add new bessel images

17 files changed, 1280 insertions, 491 deletions

diff --git a/buch/chapters/060-integral/Makefile.inc b/buch/chapters/060-integral/Makefile.inc
index 9cc5356..09be355 100644
--- a/buch/chapters/060-integral/Makefile.inc
+++ b/buch/chapters/060-integral/Makefile.inc
@@ -9,5 +9,6 @@ CHAPTERFILES = $(CHAPTERFILES)						\
 	chapters/060-integral/differentialkoerper.tex			\
 	chapters/060-integral/risch.tex					\
 	chapters/060-integral/orthogonal.tex				\
+	chapters/060-integral/legendredgl.tex				\
 	chapters/060-integral/gaussquadratur.tex			\
 	chapters/060-integral/chapter.tex
diff --git a/buch/chapters/060-integral/legendredgl.tex b/buch/chapters/060-integral/legendredgl.tex
new file mode 100644
index 0000000..9aeac40
--- /dev/null
+++ b/buch/chapters/060-integral/legendredgl.tex
@@ -0,0 +1,369 @@
+%
+% legendredgl.tex
+%
+% (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+%
+\subsection{Orthogonale Polynome und Differentialgleichungen}
+Legendre hat einen ganz anderen Zugang zu den nach ihm benannten
+Polynomen gefunden.
+Er hat sie gefunden als die Lösungen einer speziellen Differentialgleichungen.
+In diesem Abschnitt sollen diese Funktionen mit der Potenzreihen-Methode
+wiedergefunden werden.
+Dabei stellt sich heraus, dass diese Polynome auch Eigenfunktionen eines
+selbstadjungierten Differentialgoperator sind.
+Die Orthogonalität wird dann aus einer Verallgemeinerung der bekannten
+Eingeschaft folgen, dass Eigenvektoren einer symmetrischen Matrix zu 
+verschiedenen Eigenwerten orthogonal sind.
+
+\subsubsection{Legendre-Differentialgleichung}
+Die {\em Legendre-Differentialgleichung} ist die Differentialgleichung
+\begin{equation}
+(1-x^2) y'' - 2x y' + n(n+1) y = 0
+\label{buch:integral:eqn:legendre-differentialgleichung}
+\end{equation}
+für eine Funktion $y(x)$ auf dem Intervall $[-1,1]$.
+
+Sei $y(x)$ eine Lösung der Differentialgleichung
+\eqref{buch:integral:eqn:legendre-differentialgleichung}.
+Setzt man $y_s(x)=y(-x)$ in die Differentialgleichung ein, erhält
+man
+\[
+(1-x^2)y_s''(x) - 2x y'_s(x) + n(n+1)y_s(x)
+=
+(1-x^2)y''(-x) +2x y(-x) +n(n+1)y(-x).
+\]
+Ersetzt man $t=-x$, dann wird daraus
+\[
+(1-x^2)y''(t) -2t y(t) + n(n+1) y(t) = 0
+\]
+aus der Differentialgleichung
+\eqref{buch:integral:eqn:legendre-differentialgleichung}.
+Insbesondere ist die gespiegelte Funktion $y_s(x)$ ebenfalls
+eine Lösung der Differentialgleichung.
+
+Ist $y(x)$ eine Lösung der Differentialgleichung ist, dann lässt
+sie sich in die Summe einer geraden und einer ungeraden Funktion
+\[
+\left.
+\begin{aligned}
+y_g(x) &= \frac{y(x)+y(-x)}{2}\\
+y_u(x) &= \frac{y(x)-y(-x)}{2}
+\end{aligned}
+\quad
+\right\}
+\quad
+\Rightarrow
+\quad
+y(x) = y_g(x) + y_u(x)
+\]
+zerlegen, die als Linearkombinationen der beiden Lösungen
+$y(x)$ und $y_s(x)$ ebenfalls Lösungen der Differentialgleichung
+sind.
+
+\subsubsection{Potenzreihenlösung}
+Wir suchen eine Lösung in Form einer Potenzreihe um $x=0$ und 
+verwenden dazu den Ansatz
+\[
+y(x) = a_0+a_1x+a_2x^2+ \dots = \sum_{k=0}^\infty a_kx^k.
+\]
+\begin{align*}
+(1-x^2) \sum_{k=2}^\infty k(k-1)a_kx^{k-2}
+-2x\sum_{k=0}^\infty ka_kx^{k-1}
++
+n(n+1)\sum_{k=0}^\infty  a_kx^k
+&=
+0
+\\
+\sum_{k=0}^\infty (k+2)(k+1)a_{k+2}x^k
+-
+\sum_{k=2}^\infty k(k-1)a_kx^k
+-
+2\sum_{k=1}^\infty ka_kx^k
++
+n(n+1)\sum_{k=0}^\infty  a_kx^k
+&=
+0
+\end{align*}
+Die Koeffizienten zur Potenz $k$ sind daher
+\begin{align}
+k&=0:
+&
+0&=
+2a_2+n(n+1)a_0
+\notag
+\\
+&&
+a_2&=-\frac{n(n+1)}{2}a_0
+\notag
+\\
+k&=1:
+&
+0&=
+6a_3-2a_1+n(n+1)a_1
+\notag
+\\
+&&
+a_3&= \frac{2-n(n+1)}{6}a_1
+\notag
+\\
+k&>1:
+&
+0&=
+(k+2)(k+1)a_{k+2} -k(k-1)a_k -2ka_k +n(n+1) a_k
+\notag
+\\
+&&
+a_{k+2}
+&=
+\frac{ k(k+1)-n(n+1) }{(k+2)(k+1)}
+a_k
+\label{buch:integral:legendre-dgl:eqn:akrek}
+\end{align}
+Wenn $a_1=0$ und $a_0\ne 1$ ist, dann ist die Funktion $y(x)$ gerade,
+alle ungeraden Koeffizienten verschwinden.
+Ebenso verschwinden alle geraden Koeffizienten, wenn $a_0=0$ und $a_1\ne 0$.
+Für jede Lösung $y(x)$ der Differentialgleichung ist
+$y_g(x)$ ein Lösung mit $a_1=0$ und $y_u(x)$ eine Lösung mit $a_0=0$.
+Wir können die Diskussion der Lösungen daher auf gerade oder ungerade
+Lösungen einschränken.
+
+Gesucht ist jetzt eine Lösung in Form eines Polynoms.
+In diesem Fall müssen die Koeffizienten $a_k$ ab einem
+gewissen Index verschwinden.
+Dies tritt nach \eqref{buch:integral:legendre-dgl:eqn:akrek} genau
+dann auf, wenn der Zähler für ein $k$ verschwindet.
+Folglich gibt es genau dann Polynomlösungen der Differentialgleichungen,
+wenn $n$ eine natürlich Zahl ist.
+Ausserdem ist die Lösung ein Polynom $\bar{P}_n(x)$ vom Grad $n$.
+Das Polynom soll wieder so normiert sein, dass $\bar{P}_n(1)=1$ ist.
+
+Die Lösungen der Differentialgleichungen können jetzt explizit
+berechnet werden.
+Zunächst ist $\bar{P}_0(x)=1$ und $\bar{P}_1(x)=x$.
+Für $n=2$ setzen wir zunächst $a_0=1$ und $a_1=0$ und erhalten
+\[
+y(x)
+=
+1 + \frac{0(0+1) - 2(2+1)}{(0+2)(0+1)}a_0 x^2
+=
+1
+-3x^2
+\qquad\text{oder}\qquad
+\bar{P}_3(x) = \frac12(3x^2-1).
+\]
+Für $n=3$ starten wir von $a_1=1$ und $a_0=0$, was zunächst $a_2=0$
+impliziert.
+Für $a_3$ finden wir
+\[
+a_3=\frac{1(1+1)-3(3+1)}{(1+2)(1+1)} = -\frac53
+\qquad\Rightarrow\qquad
+y(x) = x-\frac53x^3
+\qquad\Rightarrow\qquad
+\bar{P}_3(x) = \frac12(5x^3-3x).
+\]
+Dies stimmt überein mit den früher gefundenen Ausdrücken für
+die Legendre-Polynome.
+
+Die Potenzreihenlösung zeigt zwar, dass es für jedes $n\in\mathbb{N}$
+eine Polynomlösung $\bar{P}_n(x)$ vom Grad $n$ gibt.
+Dies kann aber nicht erklären, warum die so gefundenen Polynome
+orthogonal sind.
+
+\subsubsection{Eigenfunktionen}
+Die Differentialgleichung
+\eqref{buch:integral:eqn:legendre-differentialgleichung}
+Kann mit dem Differentialoperator
+\[
+D = \frac{d}{dx}(1-x^2)\frac{d}{dx}
+\]
+als
+\[
+Dy + n(n+1)y = 0
+\]
+geschrieben werden.
+Tatsächlich ist
+\[
+Dy
+=
+\frac{d}{dx} (1-x^2) \frac{d}{dy}
+=
+\frac{d}{dx} (1-x^2)y'
+=
+(1-x^2)y'' -2x y'.
+\]
+Dies bedeutet, dass die Lösungen $\bar{P}_n(x)$ Eigenfunktionen
+des Operators $D$ zum Eigenwert $n(n+1)$ sind:
+\[
+D\bar{P}_n = -n(n+1) \bar{P}_n.
+\]
+
+\subsubsection{Orthogonalität von $\bar{P}_n$ als Eigenfunktionen}
+Ein Operator $A$ auf Funktionen heisst {\em selbstadjungiert}, wenn
+für zwei beliebige Funktionen $f$ und $g$ gilt
+\[
+\langle Af,g\rangle = \langle f,Ag\rangle
+\]
+gilt.
+Im vorliegenden Zusammenhang möchten wir die Eigenschaft nutzen,
+dass Eigenfunktionen eines selbstadjungierten Operatores zu verschiedenen
+Eigenwerten orthogonal sind.
+Dazu seien $Df = \lambda f$ und $Dg=\mu g$ und wir rechnen
+\begin{equation*}
+\begin{array}{rcccl}
+\langle Df,g\rangle &=& \langle \lambda f,g\rangle &=& \lambda\langle f,g\rangle
+\\
+=\langle f,Dg\rangle &=& \langle f,\mu g\rangle &=& \mu\langle f,g\rangle
+\\
+                     & &                    0   &=& (\lambda-\mu)\langle f,g\rangle
+\end{array}
+\end{equation*}
+Da $\lambda-\mu\ne 0$ ist, muss $\langle f,g\rangle=0$ sein.
+
+Der Operator $D$ ist selbstadjungiert, d.~h.
+für zwei beliebige zweimal stetig differenzierbare Funktion $f$ und $g$
+auf dem Intervall $[-1,1]$ gilt
+\begin{align*}
+\langle Df,g\rangle
+&=
+\int_{-1}^1 (Df)(x) g(x) \,dx
+\\
+&=
+\int_{-1}^1
+\biggl(\frac{d}{dx} (1-x^2)\frac{d}{dx}f(x)\biggr) g(x)
+\,dx
+\\
+&=
+\underbrace{
+\biggl[
+\biggl((1-x^2)\frac{d}{dx}f(x)\biggr) g(x)
+\biggr]_{-1}^1
+}_{\displaystyle = 0}
+-
+\int_{-1}^1
+\biggl((1-x^2)\frac{d}{dx}f(x)\biggr) \frac{d}{dx}g(x)
+\,dx
+\\
+&=
+-
+\int_{-1}^1
+\biggl(\frac{d}{dx}f(x)\biggr) \biggl((1-x^2)\frac{d}{dx}g(x)\biggr)
+\,dx
+\\
+&=
+-
+\underbrace{
+\biggl[
+f(x) \biggl((1-x^2)\frac{d}{dx}g(x)\biggr)
+\biggr]_{-1}^1}_{\displaystyle = 0}
++
+\int_{-1}^1
+f(x) \biggl(\frac{d}{dx}(1-x^2)\frac{d}{dx}g(x)\biggr)
+\,dx
+\\
+&=
+\langle f,Dg\rangle.
+\end{align*}
+Dies beweist, dass $D$ selbstadjungiert ist.
+Da $\bar{P}_n$ Eigenwerte des selbstadjungierten Operators $D$ zu
+den verschiedenen Eigenwerten $-n(n+1)$ sind, folgt auch, dass
+die $\bar{P}_n$ orthogonale Polynome vom Grad $n$ sind, die die 
+gleiche Standardierdisierungsbedingung wie die Legendre-Polyonome
+erfüllen, also ist $\bar{P}_n(x)=P_n(x)$.
+
+\subsubsection{Legendre-Funktionen zweiter Art}
+Siehe Wikipedia-Artikel \url{https://de.wikipedia.org/wiki/Legendre-Polynom}
+
+Die Potenzreihenmethode liefert natürlich auch Lösungen der
+Legendreschen Differentialgleichung, die sich nicht als Polynome
+darstellen lassen.
+Ist $n$ gerade, dann liefern die Anfangswerte $a_0=0$ und $a_1=1$ 
+eine ungerade Funktion, die Folge der Koeffizienten bricht
+aber nicht ab, vielmehr ist
+\begin{align*}
+a_{k+2}
+&=
+\frac{k(k+1)}{(k+1)(k+2)}a_k
+=
+\frac{k}{k+2}a_k.
+\end{align*}
+Durch wiederholte Anwendung dieser Rekursionsformel findet man
+\[
+a_{k}
+=
+\frac{k-2}{k}a_{k-2}
+=
+\frac{k-2}{k}\frac{k-4}{k-2}a_{k-4}
+=
+\frac{k-2}{k}\frac{k-4}{k-2}\frac{k-6}{k-4}a_{k-6}
+=
+\dots
+=
+\frac{1}{k}a_1.
+\]
+Die Lösung hat daher die Reihenentwicklung
+\[
+Q_0(x) = x+\frac13x^3 + \frac15x^5 + \frac17x^7+\dots
+=
+\frac12\log \frac{1+x}{1-x}
+=
+\operatorname{artanh}x.
+\]
+Die Funktion $Q_0(x)$ heisst {\em Legendre-Funktion zweiter Art}.
+
+Für $n=1$ wird die Reihenentwicklung $a_0=1$ und $a_1=0$ etwas
+interessanter.
+Die Rekursionsformel für die Koeffizienten ist
+\[
+a_{k+2}
+=
+\frac{k(k+1)-2}{(k+1)(k+2)} a_k.
+\qquad\text{oder}\qquad
+a_k
+=
+\frac{(k-1)(k-2)-2}{k(k-1)}
+a_{k-2}
+\]
+Man erhält der Reihe nach
+\begin{align*}
+a_2 &= \frac{-2}{2\cdot 1} a_0 = -1
+\\
+a_3 &= 0
+\\
+a_4 &= \frac{3\cdot 2-2}{4\cdot 3} a_2 = \frac{4}{4\cdot 3}a_2 = \frac13a_2 = -\frac13
+\\
+a_5 &= 0
+\\
+a_6 &= \frac{5\cdot 4-2}{6\cdot 5}a_4 = \frac{18}{6\cdot 5}a_4 = -\frac15
+\\
+a_7 &= 0
+\\
+a_8 &= \frac{7\cdot 6-2}{8\cdot 7}a_6 = \frac{40}{8\cdot 7} = -\frac17
+\\
+a_9 &= 0
+a_{10} &= \frac{9\cdot 8-2}{10\cdot 9}a_8 = \frac{70}{10\cdot 9} = -\frac19,
+\end{align*}
+woraus sich die Reihenentwicklung
+\begin{align*}
+y(x)
+&=
+-x^2 -\frac13x^4 -\frac15x^6 - \frac17x^8 -\frac19x^{10}-\dots
+\\
+&=
+-x\biggl(x+\frac13x^3 + \frac15x^5 + \frac17x^7 + \frac19x^9+\dots\biggr)
+=
+-x\operatorname{artanh}x.
+\end{align*}
+Die {\em Legendre-Funktionen zweiter Art} $Q_n(x)$  werden allerdings
+so definiert, dass gewisse Rekursionsformeln für die Legendre-Polynome,
+die wir hier nicht hergeleitet haben, auch für die $Q_n(x)$ gelten.
+In dieser Normierung muss statt des eben berechneten $y(x)$ die Funktion
+\[
+Q_1(x) = x \operatorname{artanh}x-1
+\]
+verwendet werden.
+
+
+
+
+
+
diff --git a/buch/chapters/060-integral/orthogonal.tex b/buch/chapters/060-integral/orthogonal.tex
index ceba53a..109cd61 100644
--- a/buch/chapters/060-integral/orthogonal.tex
+++ b/buch/chapters/060-integral/orthogonal.tex
@@ -508,496 +508,15 @@ 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
+%% Differentialgleichungen
 %%
-%\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}
+%\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}
+
+\input{chapters/060-integral/legendredgl.tex}
+
 \input{chapters/060-integral/gaussquadratur.tex}
diff --git a/buch/chapters/090-pde/Makefile.inc b/buch/chapters/090-pde/Makefile.inc
new file mode 100644
index 0000000..191bad6
--- /dev/null
+++ b/buch/chapters/090-pde/Makefile.inc
@@ -0,0 +1,12 @@
+#
+# Makefile.inc -- Makefile dependencies for chapter 9
+#
+# (c) 2022 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+#
+
+CHAPTERFILES = $(CHAPTERFILES)						\
+	chapters/090-pde/gleichung.tex					\
+	chapters/090-pde/separation.tex					\
+	chapters/090-pde/rechteck.tex					\
+	chapters/090-pde/kreis.tex					\
+	chapters/090-pde/chapter.tex
diff --git a/buch/chapters/090-pde/bessel/Makefile b/buch/chapters/090-pde/bessel/Makefile
new file mode 100644
index 0000000..c189517
--- /dev/null
+++ b/buch/chapters/090-pde/bessel/Makefile
@@ -0,0 +1,15 @@
+#
+# Makefile
+#
+# (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+#
+all:	besselzeros.tex besselnodes.tex pauke.pdf
+
+besselzeros.tex:	besselzeros.m
+	octave besselzeros.m
+
+besselnodes.tex:	besselnodes.m
+	octave besselnodes.m
+
+pauke.pdf:	pauke.tex besselnodes.tex
+	pdflatex pauke.tex
diff --git a/buch/chapters/090-pde/bessel/besselnodes.m b/buch/chapters/090-pde/bessel/besselnodes.m
new file mode 100644
index 0000000..0dcba3e
--- /dev/null
+++ b/buch/chapters/090-pde/bessel/besselnodes.m
@@ -0,0 +1,106 @@
+#
+# besselnodes.m
+#
+# (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+#
+global maxmu;
+maxmu = 3;
+global maxk;
+maxk = 4;
+global mu;
+
+nachkommastellen = 4;
+
+function retval = f(x)
+	global mu;
+	retval = besselj(mu, x);
+endfunction
+
+global kzeros;
+kzeros = zeros(maxk+1, maxmu+1);
+for mu = (0:maxmu)
+	k = 0;
+	x = 0.0001;
+	while (k <= maxk)
+		bracket = [ x, x+1 ];
+		if (f(bracket(1)) * f(bracket(2)) < 0) 
+			kzeros(k+1,mu+1) = fzero("f", bracket);
+			k = k + 1;
+		endif
+		x = x + 1;
+	endwhile
+endfor
+
+xshift = 4;
+yshift = 4;
+global	r;
+r = 1.8;
+
+function	retval = anderefarbe(f)
+	if (1 == strcmp("red", f))
+		retval = "blue";
+	else
+		retval = "red";
+	endif
+endfunction
+
+function sektor(fn, mu, k, w0, w1, startfarbe)
+	global	kzeros;
+	global 	r;
+	fprintf(fn, "\\begin{scope}\n");
+	fprintf(fn, "\\clip (0,0)--(%.4f:%.4f) arc (%.4f:%.4f:%.4f)--cycle;\n",
+		w0, r, w0, w1, r);
+	faktor = kzeros(k+1,mu+1);
+
+	K = k + 1;
+	farbe = startfarbe;
+	while (K > 0)
+		R = r * kzeros(K, mu+1) / faktor;
+		fprintf(fn, "\\fill[color=%s!20] ", farbe);
+		fprintf(fn, "(0,0) circle[radius=%.4f];\n", R);
+		farbe = anderefarbe(farbe);
+		K = K-1;
+	end
+	fprintf(fn, "\\end{scope}\n");
+endfunction
+
+fn = fopen("besselnodes.tex", "w");
+
+#fprintf(fn, "\\begin{tikzpicture}[>=latex,thick]\n");
+
+for mu = (0:maxmu)
+	if (mu > 0)
+		winkel = 180 / mu;
+	else
+		winkel = 360;
+	endif
+	for k = (0:maxk)
+		fprintf(fn, "\\begin{scope}[xshift=%.3fcm,yshift=-%.3fcm]\n",
+			mu * xshift, k * yshift);
+		for w0 = (0:2*winkel:360)
+			sektor(fn, mu, k, w0, w0 + winkel, "red");
+			if (winkel < 270)
+				sektor(fn, mu, k, w0 + winkel, w0 + 2 * winkel, "blue");
+			endif
+		endfor
+
+		fprintf(fn, "\\draw (0,0) circle[radius=%.4f];\n", r);
+
+		fprintf(fn, "\\end{scope}\n\n");
+	endfor
+endfor
+
+for mu = (0:maxmu)
+	fprintf(fn, "\\node at (%.4f,%.4f) [above] {$\\mu=%d$};\n",
+		mu * xshift, 0.5 * yshift, mu);
+endfor
+
+for k = (0:maxk)
+	fprintf(fn, "\\node at (%.4f,%.4f) [above,rotate=90] {$k=%d$};\n",
+		-0.5 * xshift, -k * yshift, k);
+endfor
+
+#fprintf(fn, "\\end{tikzpicture}\n");
+
+fclose(fn);
+
diff --git a/buch/chapters/090-pde/bessel/besselzeros.m b/buch/chapters/090-pde/bessel/besselzeros.m
new file mode 100644
index 0000000..9c8fa9d
--- /dev/null
+++ b/buch/chapters/090-pde/bessel/besselzeros.m
@@ -0,0 +1,70 @@
+#
+# besselzeros.m -- find zeros of bessel functions
+#
+# (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+# 
+global maxmu;
+maxmu = 7;
+global maxk;
+maxk = 10;
+global mu;
+
+nachkommastellen = 4;
+
+function retval = f(x)
+	global mu;
+	retval = besselj(mu, x);
+endfunction
+
+kzeros = zeros(maxk+1, maxmu+1);
+for mu = (0:maxmu)
+	k = 0;
+	if (mu > 0)
+		kzeros(1, mu+1) = 0;
+		k = k+1;
+	endif
+	x = 0.0001;
+	while (k <= maxk)
+		bracket = [ x, x+1 ];
+		if (f(bracket(1)) * f(bracket(2)) < 0) 
+			kzeros(k+1,mu+1) = fzero("f", bracket);
+			k = k + 1;
+		endif
+		x = x + 1;
+	endwhile
+endfor
+
+# kzeros
+
+fn = fopen("besselzeros.tex", "w");
+
+fprintf(fn, "\\begin{tabular}{|>{$}c<{$}");
+for mu = (0:maxmu)
+	fprintf(fn, "|>{$}r<{$}");
+endfor
+fprintf(fn, "|}\n");
+
+fprintf(fn, "\\hline\n");
+fprintf(fn, " k ");
+for mu = (0:maxmu)
+	fprintf(fn, "& \\mu = %d ", mu);
+endfor
+fprintf(fn, "\\\\\n");
+fprintf(fn, "\\hline\n");
+
+for k = (0:maxk)
+	fprintf(fn, " %d ", k);
+	for mu = (0:maxmu)
+		value = kzeros(k+1, mu+1);
+		if (value == 0) 
+			fprintf(fn, "& 0\\phantom{.%0*d}", nachkommastellen, 0);
+		else
+			fprintf(fn, "& %*.*f", nachkommastellen+4, nachkommastellen, kzeros(k+1, mu+1));
+		endif
+	endfor
+	fprintf(fn, "\\\\\n");
+endfor
+fprintf(fn, "\\hline\n");
+fprintf(fn, "\\end{tabular}\n");
+
+fclose(fn);
diff --git a/buch/chapters/090-pde/bessel/besselzeros.tex b/buch/chapters/090-pde/bessel/besselzeros.tex
new file mode 100644
index 0000000..1b8d33b
--- /dev/null
+++ b/buch/chapters/090-pde/bessel/besselzeros.tex
@@ -0,0 +1,17 @@
+\begin{tabular}{|>{$}c<{$}|>{$}r<{$}|>{$}r<{$}|>{$}r<{$}|>{$}r<{$}|>{$}r<{$}|>{$}r<{$}|>{$}r<{$}|>{$}r<{$}|}
+\hline
+ k & \mu = 0 & \mu = 1 & \mu = 2 & \mu = 3 & \mu = 4 & \mu = 5 & \mu = 6 & \mu = 7 \\
+\hline
+ 0 &   2.4048& 0\phantom{.0000}& 0\phantom{.0000}& 0\phantom{.0000}& 0\phantom{.0000}& 0\phantom{.0000}& 0\phantom{.0000}& 0\phantom{.0000}\\
+ 1 &   5.5201&   3.8317&   5.1356&   6.3802&   7.5883&   8.7715&   9.9361&  11.0864\\
+ 2 &   8.6537&   7.0156&   8.4172&   9.7610&  11.0647&  12.3386&  13.5893&  14.8213\\
+ 3 &  11.7915&  10.1735&  11.6198&  13.0152&  14.3725&  15.7002&  17.0038&  18.2876\\
+ 4 &  14.9309&  13.3237&  14.7960&  16.2235&  17.6160&  18.9801&  20.3208&  21.6415\\
+ 5 &  18.0711&  16.4706&  17.9598&  19.4094&  20.8269&  22.2178&  23.5861&  24.9349\\
+ 6 &  21.2116&  19.6159&  21.1170&  22.5827&  24.0190&  25.4303&  26.8202&  28.1912\\
+ 7 &  24.3525&  22.7601&  24.2701&  25.7482&  27.1991&  28.6266&  30.0337&  31.4228\\
+ 8 &  27.4935&  25.9037&  27.4206&  28.9084&  30.3710&  31.8117&  33.2330&  34.6371\\
+ 9 &  30.6346&  29.0468&  30.5692&  32.0649&  33.5371&  34.9888&  36.4220&  37.8387\\
+ 10 &  33.7758&  32.1897&  33.7165&  35.2187&  36.6990&  38.1599&  39.6032&  41.0308\\
+\hline
+\end{tabular}
diff --git a/buch/chapters/090-pde/bessel/pauke.pdf b/buch/chapters/090-pde/bessel/pauke.pdf
new file mode 100644
index 0000000..54edc20
--- /dev/null
+++ b/buch/chapters/090-pde/bessel/pauke.pdf
diff --git a/buch/chapters/090-pde/bessel/pauke.tex b/buch/chapters/090-pde/bessel/pauke.tex
new file mode 100644
index 0000000..bba092e
--- /dev/null
+++ b/buch/chapters/090-pde/bessel/pauke.tex
@@ -0,0 +1,21 @@
+%
+% pauke.tex -- template for standalon tikz images
+%
+% (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}
+\begin{tikzpicture}[>=latex,thick,scale=\skala]
+
+\input{besselnodes.tex}
+
+\end{tikzpicture}
+\end{document}
+
diff --git a/buch/chapters/090-pde/chapter.tex b/buch/chapters/090-pde/chapter.tex
new file mode 100644
index 0000000..543a92d
--- /dev/null
+++ b/buch/chapters/090-pde/chapter.tex
@@ -0,0 +1,30 @@
+%
+% chapter.tex -- Kapitel zu partiellen Differentialgleichungen
+%
+% (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+%
+% !TeX spellcheck = de_CH
+\chapter{Partielle Differentialgleichungen
+\label{buch:chapter:pde}}
+\lhead{Partielle Differentialgleichungen}
+\rhead{}
+Partielle Differentialgleichungen sind eine besonders ergiebige
+Quelle für Anwendungen spezieller Funktionen.
+Die Separationsmethode zum Beispiel für die Wellengleichung
+auf gewissen, besonders einfachen Gebieten wie Rechtecken,
+Kreisscheiben oder Kugel führt auf gewöhnliche Differentialgleichungen,
+deren Lösungen spezielle Funktionen sind.
+
+\input{chapters/090-pde/gleichung.tex}
+\input{chapters/090-pde/separation.tex}
+\input{chapters/090-pde/rechteck.tex}
+\input{chapters/090-pde/kreis.tex}
+
+%\section*{Übungsaufgaben}
+%\rhead{Übungsaufgaben}
+%\aufgabetoplevel{chapters/020-exponential/uebungsaufgaben}
+%\begin{uebungsaufgaben}
+%\uebungsaufgabe{0}
+%\uebungsaufgabe{1}
+%\end{uebungsaufgaben}
+
diff --git a/buch/chapters/090-pde/gleichung.tex b/buch/chapters/090-pde/gleichung.tex
new file mode 100644
index 0000000..07dd2ff
--- /dev/null
+++ b/buch/chapters/090-pde/gleichung.tex
@@ -0,0 +1,11 @@
+%
+% gleichung.tex
+%
+% (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+%
+\section{Gleichungen und Randbedingungen
+\label{buch:pde:section:gleichungen-und-randbedingungen}}
+
+\subsection{Laplace-Operator}
+
+\subsection{Orthogonalität}
diff --git a/buch/chapters/090-pde/kreis.tex b/buch/chapters/090-pde/kreis.tex
new file mode 100644
index 0000000..a54ce38
--- /dev/null
+++ b/buch/chapters/090-pde/kreis.tex
@@ -0,0 +1,219 @@
+%
+% kreis.tex
+%
+% (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+%
+\section{Kreisförmige Membran
+\label{buch:pde:section:kreis}}
+In diesem Abschnitt soll die Differentialgleichung einer kreisförmigen
+Membran mit Hilfe der Separationsmethode gelöst werden.
+Dabei werden die Bessel-Funktionen als Lösungsfunktionen 
+auftreten und die Eigenfrequenzen werden durch ihre Nullstellen
+berechnet.
+
+\subsection{Differentialgleichung und Randbedingung}
+Die Wellengleichung auf einem Kreisgebiet mit Radius $r_0$ 
+lässt sich am besten mit Hilfe von Polarkoordinaten $(r,\varphi)$
+ausdrücken.
+Gesucht ist also eine Funktion $u(t,r,\varphi)$ gesucht, wobei
+$0\le r<r_0$ und $0\le \varphi\le 2\pi$.
+Die Funktion muss eine Lösung der Wellengleichung
+\[
+\frac{1}{c^2}\frac{\partial^2u}{\partial t^2} = \Delta u
+\]
+sein.
+
+Der Laplace-Operator hat in Polarkoordinaten die Form
+\begin{equation}
+\Delta
+=
+\frac{\partial^2}{\partial r^2}
++
+\frac1r
+\frac{\partial}{\partial r}
++
+\frac{1}{r 2}
+\frac{\partial^2}{\partial\varphi^2}.
+\label{buch:pde:kreis:laplace}
+\end{equation}
+Die Differentialgleichung ist 
+\[
+\frac{1}{c^2} \frac{\partial^2 u}{\partial t^2}
+=
+\Delta u.
+\]
+Die Separation der Zeit führt auf die Eigenwertgleichung
+\[
+\Delta U(r,\varphi) = -\lambda^2 U(r,\varphi)
+\]
+für eine Funktion, die nur von $r$ und $\varphi$ abhängt.
+
+Die Randbedingungen besagen, dass $u(t,r_0,\varphi)=0$ für $t>0$.
+Dies bedeutet für die Funktion $U(r,\varphi)$, dass
+$U(r_0,\varphi)=0$ sein muss für alle $\varphi$.
+
+Die Bedingungen an $U$ reichen aber nicht ganz.
+Alle Koordinaten $(0,\varphi)$ bezeichnen ja gleichermassen
+den Nullpunkt des Koordinatensystems, es muss also auch sichergestellt
+sein, dass $U(0,\varphi)$ für alle $\varphi$ den gleichen Wert gibt.
+
+\subsection{Separation}
+Das Eigenwertproblem $\Delta U=-\lambda^2 U$ soll jetzt in Polarkoordinaten
+separiert werden.
+Dazu schreiben wir die Lösung als
+\[
+U(r,\varphi)
+=
+R(r)\cdot \Phi(\varphi).
+\]
+Die Randbedingungen an $U$ werden zu $R(r_0)=0$.
+
+Im Ursprung des Koordinatensystems ist die Randbedingung etwas
+komplizierter.
+Wenn $R(0)=0$ ist, dann ist sichergestellt, dass
+$U(0,\varphi)=R(0)\Phi(\varphi)0$ ist, dass also der Wert unabhängig
+ist von $\varphi$.
+Wenn aber $R(0)\ne 0$ ist, dann kann die geforderte Unabhängigkeit
+von $\varphi$ nur erfüllt werden, wenn $\Phi(\varphi)$ konstant ist.
+Da die Funktion aber auch noch differenzierbar sein soll, darf es
+an der Stelle $r=0$ keine ``Spitze'' geben, die Ableitung $R'(0)$
+muss also auch $=0$ sein.
+% XXX Evtl Bild zur Illustration dieses Problems
+
+Die Differntialgleichungen wird mit der Form~\eqref{buch:pde:kreis:laplace}
+des Laplace-Operators
+\[
+\Delta U
+=
+R''(r) \Phi(\varphi)
++
+\frac1r R'(r)\Phi(\varphi)
++
+\frac{1}{r^2} R(r)\Phi''(\varphi)
+=
+-\lambda^2
+R(r)\Phi(\varphi)
+\]
+Nach Division durch die rechte Seite erhalten wir
+\[
+\frac{R''(r)}{R(r)}
++
+\frac1r \frac{R'(r)}{R(r)}
++
+\frac{1}{r^2} \frac{\Phi''(\varphi)}{\Phi(\varphi)}
+=
+-\lambda^2
+\]
+Im letzten Term auf der linken Seite kommen die Variablen $r$ und $\varphi$
+gemischt vor, man muss also die Gleichung erst mit $r^2$ multiplizieren,
+bevor man sie in 
+\[
+\frac{r^2R''(r)+rR'(r)+\lambda^2 r^2R(r)}{R(r)}
+=
+-\frac{\Phi''(\varphi)}{\Phi(\varphi)}
+\]
+separieren kann.
+Die beiden Seiten sind also konstant, wir nennen die gemeinsame
+Konstanten $\mu^2$, das vereinfacht die Lösung der Gleichung
+für $\Phi(\varphi)$.
+
+Die Gleichung für $\Phi$ hat für $\mu\ne 0$ die Lösungen
+\begin{align*}
+\Phi(\varphi) &= \cos\mu\varphi
+\text{und}\qquad
+\Phi(\varphi) &= \sin\mu\varphi.
+\end{align*}
+Die Lösung muss aber auch stetig sein, d.~h.~es muss $\Phi(0)=\Phi(2\pi)$
+gelten.
+Dies ist nur möglich, wenn $\mu$ eine ganze Zahl ist.
+
+Für $\mu=0$ hat das charakteristische Polynome eine doppelte Nullstelle,
+die allgemeine Lösung lautet daher
+\[
+\Phi(\varphi)= C \varphi + D.
+\]
+Die Funktion $\Phi$ muss aber auch stetig sein, d.~h.~$\Phi(0)=\Phi(2\pi)$,
+das ist mit $C\ne 0$ nicht möglich, somit kommt für $\mu=0$ nur die
+Lösung $\Phi(\varphi)=D$ in Frage.
+
+Die Gleichung für $R(r)$ wird jetzt
+\begin{equation}
+r^2R''(r) + rR'(r)+(\lambda^2 r^2-\mu^2)R(r)
+=
+0.
+\label{buch:pde:kreis:Rdgl}
+\end{equation}
+Bis auf den Faktor $\lambda^2$ ist dies eine Besselsche Differentialgleichung.
+
+\subsection{Umformung in eine Besselsche Differentialgleichung}
+Die Funktion $y(x) = J_\mu(sx)$ hat die Ableitungen
+\begin{align*}
+y'(x) &= sJ'_mu(sx)      
+\\
+y''(x) &= s^2J''_\mu(sx)
+\end{align*}
+Setzt man dies in die Besselsche Differentialgleichung für $J_\mu$ an
+der Stelle $sx$ ein, erhält man
+\[
+s^2x^2 J''_\mu(sx) + sx J'_\mu(sx) + (s^2x^2 -\mu^2) J_\mu(sx) = 0.
+\]
+Die Differentialgleichung \eqref{buch:pde:kreis:Rdgl} der Funktion $R(r)$
+wird also gelöst von den Funktionen $R(r) = J_\mu(\lambda r)$.
+
+\subsection{Eigenfrequenzen}
+Im vorangegangenen Abschnitt haben wir gefunden, dass die Lösungen
+für $R(r)$ die Funktionen $J_\mu(\lambda r)$ sind.
+Bis jetzt haben wir aber nicht nachgeprüft, dass die Randbedingung
+eingehalten wird. 
+Diese ist erfüllt, $R(r_0)=0$ ist.
+Es muss also
+$J_\mu(\lambda r_0)=0$ sein, oder $\lambda r_0$ muss eine
+Nullstelle von $J_{\mu}$ sein.
+Bezeichnen wir die Nullstellen von $J_\mu$ mit $j_{\mu k}$, wobei $k$
+eine natürliche Zahl ist, dann muss
+\[
+\lambda = \frac{j_{\mu k}}{r_0}
+\]
+sein.
+Die Eigenfrequenzen der kreisförmigen Membran werden also im Wesentlichen
+durch die Nullstellen der Bessel-Funktionen gegeben.
+
+Zu jedem ganzzahligen $\mu$ gibt es also eine Folge $j_{\mu k}/r_0$  von
+Eigenfrequenzen. 
+Die Lösungen mit Index $k$ der Differentialgleichung mit Index $k$ hat
+die Form
+\[
+U_{\mu k}(r,\varphi)
+=
+C \cos(\mu \varphi+\delta)
+J_{\mu}\biggl(
+\frac{j_{\mu k}}{r_0}r
+\biggr)
+\]
+Der Faktor $J_{\mu}$ hat $k$ weitere Nullstellen für Radien $r<r_0$,diese
+gehören zu kreisförmigen Knotenlinien der Membran, dort bewegt sie sich
+nicht.
+Der Faktor $\cos(\mu\varphi+\delta)$ hat $2\mu$ Nullstellen im Intervall
+$[0,2\pi)$, es gibt also noch zusätzlich $\mu$ diametrale Knotenlinien.
+Nur für $\mu=0$ gibt es Lösungen, die keine radialen Knotenlinien haben,
+da in diesem Fall $\Phi$ eine konstante Funktion sein muss.
+
+\begin{table}
+\centering
+\input{chapters/090-pde/bessel/besselzeros.tex}
+\caption{Nullstellen der Bessel-Funktionen
+\label{buch:pde:kreis:table:besselzeros}}
+\end{table}
+
+\begin{figure}
+\centering
+\includegraphics[width=\textwidth]{chapters/090-pde/bessel/pauke.pdf}
+\caption{Vorzeichen der Lösungsfunktionen und Knotenlinien
+für verschiedene Werte von $\mu$ und $k$.
+Die Bereiche, in denen die Lösungsfunktion positiv sind, ist 
+rot dargestellt, die negativen Bereiche blau.
+In jeder Darstellung gibt es genau $k+\mu$ Knotenlinien.
+Die Radien der kreisförmigen Knotenlinien müssen aus den Nullstellen
+der Besselfunktionen berechnet werden.
+\label{buch:pde:kreis:fig:pauke}}
+\end{figure}
diff --git a/buch/chapters/090-pde/rechteck.tex b/buch/chapters/090-pde/rechteck.tex
new file mode 100644
index 0000000..944fbf1
--- /dev/null
+++ b/buch/chapters/090-pde/rechteck.tex
@@ -0,0 +1,7 @@
+%
+% rechteck.tex
+%
+% (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+%
+\section{Rechteckige Membran
+\label{buch:pde:section:rechteck}}
diff --git a/buch/chapters/090-pde/separation.tex b/buch/chapters/090-pde/separation.tex
new file mode 100644
index 0000000..81195d3
--- /dev/null
+++ b/buch/chapters/090-pde/separation.tex
@@ -0,0 +1,391 @@
+%
+% separation.tex
+%
+% (c) 2021 Prof Dr Andreas Müller, OST Ostschweizer Fachhochschule
+%
+\section{Separationsmethode
+\label{buch:pde:section:separation}}
+Die Existenz der Lösung einer gewöhnlichen Differentialgleichung
+ist unter einigermassen milden Bedingungen in der Nähe der
+Anfangsbedingung garantiert.
+Ausserdem steht eine ganze Reihe von Lösungsverfahren zur
+Verfügung, nicht zuletzt das Potenzreihenverfahren, welches in
+Kapitel~\ref{buch:chapter:differential} beschrieben wurde.
+Das Ziel dieses Abschnitts ist eine Methode vorzustellen, mit
+der partielle Differentialgleichungen auf gewöhnliche
+Differentialgleichungen zurückgeführt werden können.
+
+%
+% Ansatz
+%
+\subsection{Separationsansatz}
+Die Separationsmethode ist motiviert durch die Beobachtung, dass in
+vielen partiellen Differentialgleichungen die Ableitungen nach
+verschiedenen Variablen sich in verschiedenen Termen befinden und
+sich daher algebraisch trennen lassen.
+Für eine beliebige Funktion bringt das nicht viel, aber für
+Funktionen mit einer speziellen Form kann man daraus eine Vereinfachung
+ableiten.
+
+%
+% Prinzip der Separation
+%
+\subsubsection{Prinzip}
+Die Grundlage der Separationsmethode ist die Idee, die Differentialgleichung
+in zwei Teile aufzuteilen, die keine gemeinsamen Variablen enthalten.
+Eine partielle Differentialgleichungen in einem zweidimensionalen
+Gebiet mit den Koordinaten $x$ und $y$ soll so umgeformt
+werden, dass auf der linken Seite des Gleichheitszeichens nur
+die Variable $x$ vorkommt und auf der rechten nur die Variable $y$.
+Es entsteht also eine Gleichung der Form
+\begin{equation}
+F(x) = G(y).
+\label{buch:pde:ansatz:eqn:F=G}
+\end{equation}
+Wie so etwas gehen gehen kann wird weiter unten untersucht.
+
+Betrachtet hält man in der Gleichung~\eqref{buch:pde:ansatz:eqn:F=G}
+die Variable $x$ fest, steht links eine fest Zahl, schreiben wir 
+sie $\lambda$.
+Die Gleichung wird also zu
+\[
+\lambda = G(y),
+\]
+sie muss für alle $y$ gelten.
+Es folgt dann, dass die rechte Seite gar nicht von $y$ abhängen kann.
+Für jeden Wert von $y$ muss $G$ den gleichen Wert $\lambda$ geben.
+
+Wenn aber $G$ konstant ist und immer den Wert $\lambda$ ergibt, dann
+ist die Gleichung~\eqref{buch:pde:ansatz:eqn:F=G} auch gleichbedeutend
+mit der Gleichung
+\[
+F(x) = \lambda,
+\]
+$F$ muss also auch konstant sein.
+
+Die algebraische Trennung der beiden Variablen $x$ und $y$ hat also 
+zur Folge, dass die beiden Seiten der Gleichung gar nicht varieren
+können, beide Seiten müssen konstant sein.
+Die Konstante ist allerdings nicht bekannt und muss im Laufe der
+weiteren Lösungsschritte der Gleichung bestimmt werden.
+
+Die Überlegungen funktionieren auch für eine grössere Zahl von
+Variablen.
+Entscheidend ist nur, dass die einen Variablen, zum Beispiel
+$x_1,\dots,x_k$, nur auf der linken Seite vorkommen und die anderen,
+wir nennen sie $x_{k+1},\dots,x_n$ nur auf der rechten.
+Die Gleichung hat dann die Form
+\begin{equation}
+F(x_1,\dots,x_k)
+=
+G(x_{k+1},\dots,x_n).
+\label{buch:pde:ansatz:eqn:FF=GG}
+\end{equation}
+Setzt man feste Werte von $x_1,\dots,x_k$ ein, ist die linke Seite
+eine Zahl, die wir wieder $\lambda$ nennen können.
+Es muss also für alle $x_{k+1},\dots,x_n$ gelten, dass
+$G(x_{k+1},\dots,x_n)=\lambda$ ist.
+Daher ist $G$ eine Konstante, sie ist gar nicht von den Variablen
+abhängig.
+Wenn aber die rechte Seite konstant ist, dann muss auch für alle
+$x_1,\dots,x_k$ gelten, dass $F(x_1,\dots,x_k)=\lambda$ ist,
+die linke Seite kann also auch nicht varieren.
+
+\begin{prinzip}
+In einer Gleichung
+\[
+F(x_1,\dots,x_k) = G(x_{k+1},\dots,x_n),
+\]
+in der die linke Seite nur von $x_1,\dots,x_k$ abhängt und die
+rechte nur von $x_{k+1},\dots,x_n$ müssen beide Seiten konstant sein.
+\end{prinzip}
+
+%
+% Beispiel zur Erklärung des Separationsvorgehens
+%
+\subsubsection{Ein Beispiel}
+In der Differentialgleichung
+\[
+x\frac{\partial u}{\partial x}
+-
+y^2\frac{\partial^2 u}{\partial y^2}
+=
+y^4
+\]
+kommen die Ableitungen nach $x$ und $y$ in verschiedenen Termen vor.
+Wir versuchen daher, auch die Lösungsfunktion als Summe
+\[
+u(x,y) = X(x) + Y(y)
+\]
+von Termen zu schreiben, die nur von jeweils einer Variablen abhängen.
+Setzt man dies in die Differentialgleichung ein, erhält man
+\[
+x\frac{\partial}{\partial x}(X(x)+Y(y))
+-y^2\frac{\partial}{\partial y}(X(x)+Y(y))
+=
+xX'(x) -y^2Y'(y)
+=
+y^4.
+\]
+Indem man den Term mit $y$ auf die rechte Seite schafft, findet man
+die Gleichung
+\[
+xX'(x) = y^2Y'(y) + y^4,
+\]
+in der die Variablen $x$ und $y$ separiert sind.
+Es folgt, dass beide Seiten konstant sein müssen, es gibt also eine
+Konstante $\lambda$ derart, dass
+\[
+xX'(x) = \lambda
+\qquad\text{und}\qquad
+y^2Y''(y) +y^4 = \lambda.
+\]
+Diese beiden Gleichungen lassen sich als Differentialgleichungen in
+der üblicheren Form als
+\begin{align*}
+X'(x) &= \frac{\lambda}{x}
+&&\Rightarrow&
+X(x) &= \int \frac{\lambda}{x}\,dx = \lambda \log x + C
+\\
+Y''(y) &= \frac{\lambda - y^4}{y^2}
+&&\Rightarrow&
+Y'(y)
+&=
+\int \frac{\lambda-y^4}{y^2}\,dy
+=
+-\frac{\lambda}{y}-\frac{y^3}3 + D
+\\
+&
+&&\Rightarrow&
+Y(y)
+&=
+\int Y'(y)\,dy
+=
+-\lambda \log y - \frac{y^4}{12} +Dy +E
+\end{align*}
+schreiben und im Falle von $X(x)$ mit einem Integral lösen.
+$Y(y)$ benötigt zwei Integrationen, ist aber ansonsten nicht
+schwieriger zu bestimmen.
+
+Das Beispiel zeigt, dass ein Separationsansatz ermöglicht, eine
+partielle Differntialgleichung in mehrere gewöhnliche Differentialgleichungen
+zu zerlegen, eine für jede Variable, und zu lösen.
+
+%
+% Anpassung des Ansatzes an die Randbedingungen
+%
+\subsubsection{Separationsansatz und Randbedingungen}
+Die im Beispiel gewählte Aufteilung der Lösungsfunktion in eine
+Summe macht es sehr schwierig, Randbedingungen der partiellen
+Differentialgleichungen in Randbedingungen der gewöhnlichen
+Differentialgleichungen zu übersetzen.
+
+Als Beispiel dieser Schwierigkeit betrachten wir die Differentialgleichung
+\[
+\Delta u
+=
+\frac{\partial^2 u}{\partial x^2}
++
+\frac{\partial^2 u}{\partial y^2}
+=
+a u
+\]
+auf dem Gebiet
+$\Omega = [0,a]\times [0,b] = \{(x,y)\in\mathbb{R}^2\mid 0<x<a\wedge 0<y<b\}$
+mit den Randwerten $u(x,y)=0$ für Punkte auf dem Rand von $\Omega$.
+Genauer:
+\[
+\begin{aligned}
+u(0,y) &= 0,& u(a,y) &= 0&&\text{für $0<y<b$} \\
+u(x,0) &= 0,& u(x,b) &= 0&&\text{für $0<x<a$}.
+\end{aligned}
+\]
+Ein Ansatz der Form $u(x,y)=X(x) + Y(y)$ bedeutet für die
+Randwerte $u(x,y)=0$, dass auf dem Rand $X(x)=-Y(y)$ gelten muss.
+Das bedeutet aber, dass $X(0) = -Y(y)$, $Y$ müsste also konstant
+sein.
+
+Ein Produktansatz löst das Problem.
+Wir verwenden stattdessen einen Produktansatz
+$u(x,y) = X(x)\cdot Y(y)$, wobei die Funktionen $X(x)$ und $Y(y)$
+nicht konstant sein sollen.
+Die Randbedingungen sind
+\[
+\begin{aligned}
+u(0,y) &= X(0) Y(y) = 0&&\Rightarrow& X(0)&=0\\
+u(a,y) &= X(a) Y(y) = 0&&\Rightarrow& X(a)&=0\\
+u(x,0) &= X(x) Y(0) = 0&&\Rightarrow& Y(0)&=0\\
+u(x,b) &= X(x) Y(b) = 0&&\Rightarrow& Y(b)&=0.
+\end{aligned}
+\]
+Der Produktansatz ermöglicht also, die Randbedingungen für die Funktion
+$u(x,y)$ in Randbedingungen für die Funktionen $X(x)$ oder $Y(y)$
+umzuwandeln.
+
+%
+% Eigenwertprobleme
+%
+\subsection{Eigenwertproblem}
+Viele partielle Differentialgleichungen der mathematischen Physik
+sind zeitabhängig, aber das räumliche Gebiet, in dem sie 
+definiert sind, ist nicht von der Zeit abhängig.
+Dies 
+
+\subsubsection{Wellengleichung}
+Die Schwingung einer ebenen Membran, die in ein emGebiet
+$G\subset\mathbb{R}^n$ eingespannt ist, wird durch die
+Wellengleichung
+\begin{equation}
+\frac{1}{c^2} \frac{\partial^2 u}{\partial t^2} = \Delta u,
+\label{buch:pde:separation:wellengleichung}
+\end{equation}
+beschrieben.
+Darin ist $u(t,x)$ die Auslenkung der Membran zur Zeit $t>0$ in einem
+Punkt $x\in G$ des Gebietes $G$ ist.
+Die Randbedingungen zerfallen in zwei Teile:
+\begin{itemize}
+\item
+Bedingungen, die wiedergeben, dass die Membran in einen 
+Rahmen eingespannt und damit unbeweglich ist.
+Dies bedeutet, dass $u(t,x)=0$ für alle Zeiten $t>0$ und für 
+Randpunkte $x\in\partial G$ von $G$ ist.
+\item
+Bedingungen, die Auslenkung und Geschwindigkeit der Membran zur
+Zeit $t=0$ beschreiben, typischerweise ind er Form
+\begin{align*}
+u(0,x) = f(x),
+\frac{\partial u}{\partial t}(0,x) = g(x)
+\end{align*}
+wobei $f(x)$ und $g(x)$ Funktionen auf dem Gebiet $G$ sind.
+\end{itemize}
+
+In der Zeitableitung auf der linken Seite
+von~\eqref{buch:pde:separation:wellengleichung}
+kommen die Ortskoordinaten nicht vor und im Laplace-Operator
+auf der rechten Seite tritt die Zeit nicht auf.
+Es ist daher naheliegend zu versuchen, die Lösung der Differntialgleichung
+als Produkt
+\[
+u(t,x) = T(t) \cdot U(x)
+\]
+zu schreiben.
+Wendet man die Differentialgleichung darauf an, wird daraus die Gleichung
+\[
+\frac{1}{c^2}
+T''(t)\cdot U(x)
+=
+T(t) \cdot \Delta U(x).
+\]
+Indem man druch $T(t)$ und $U(x)$ teilt, entsteht die separierte Gleichung
+\[
+\frac{1}{c^2} \frac{T''(t)}{T(t)}
+=
+\frac{\Delta U(x)}{U(x)}.
+\]
+Die linke Seite ist nur von der Zeit abhängig, die rechte nur von den
+Ortskoordinaten.
+Damit ist die Differentialgleichung separiert und das Problem darauf
+reduziert, die gewöhnliche Differentialgleichung 
+\[
+T''(t) = \lambda T(t)
+\]
+und die partielle Differentialgleichung
+\[
+\Delta U(x) = \lambda U(x)
+\]
+niedrigerer Dimension zu lösen.
+
+\subsubsection{Allgemeine Situation}
+Das Definitionsgebiet der partiellen Differentialgleichung ist 
+also von der Form $\mathbb{R}^+\times G$, wobei $G\subset\mathbb{R}^n$
+ein räumliches Gebiet ist und $\mathbb{R}^+$ die Zeitachse.
+Auch die Randbedingungen zerfallen in zwei Arten:
+\begin{itemize}
+\item
+Bedingungen über die Lösungsfunktion zur Zeit $t=0$ im inneren des
+räumliche Gebietes $G$, zum Beispiel
+die Anfangsauslenkung und/oder Anfangsgeschwindigkeit einer schwingenden
+Saite oder Membran.
+\item
+Bedingungen über die Lösungsfunktion auf dem Rand $\partial G$ von
+$G$ für alle Zeiten $t>0$, zum Beispiel die Bedingung, dass die
+Membran fest eingespannt ist.
+\end{itemize}
+Oft zerfällt auch der Differentialoperator in Zeitableitungen
+und einen zeitunabhängigen Teil der nur Ableitungen nach den
+Ortsvariablen enthält.
+Die Wellengleichung
+\[
+\frac{1}{c^2}
+\frac{\partial^2}{\partial t^2} u
+=
+\Delta u
+\qquad\Leftrightarrow\qquad
+\biggl(
+\frac{1}{c^2}\frac{\partial^2}{\partial t^2} - \Delta
+\biggr) u = 0
+\]
+enthält Ableitungen nach der Zeit, die nicht von Ortskoordinaten
+abhängig sind.
+Der Laplace-Operator $\Delta$ ist nicht von der Zeitabhängig und das
+Gebiet $G$ hängt ebenfalls nicht von der Zeit ab.
+
+\subsubsection{Separation der Zeit}
+Unter den gegeben Voraussetzungen ist es naheliegend, die Lösungsfunktion
+$u(t,x)$ als Produkt
+\[
+u(t,x) = T(t) \cdot U(x),\qquad t\in\mathbb{R}^+, x\in G
+\]
+anzusetezen.
+Die Wellengleichung wird dann
+\[
+\frac{1}{c^2}
+T''(t)\cdot U(x)
+=
+T(t)\cdot\Delta U(x)
+\]
+und nach Separation
+\[
+\frac{1}{c^2} \frac{T''(t)}{T(t)}
+=
+\frac{\Delta U(x)}{U(x)}.
+\]
+Es gibt also eine gemeinsame Konstante.
+Da wir Schwingungslösungen erwarten, für die $T''(t) = -\omega^2 T(t)$
+ist, schreiben wir die gemeinsame Konstante als $-\lambda^2$, was
+später die Formeln vereinfachen wird.
+Die separierten Differentialgleichungen werden jetzt
+\begin{align*}
+\frac{1}{c^2}
+\frac{T''(t)}{T(t)}
+&=
+-\lambda^2
+&&\Rightarrow&
+T''(t)-c^2\lambda T(t)&=0
+&&\Rightarrow&
+T''(t) &= A \cos(c\sqrt\lambda t) + B \sin(c \lambda t)
+\\
+&&&&&&&&
+       &= C \cos(c \lambda t+\delta)
+\\
+\frac{\Delta U(x)}{U(x)}&=-\lambda^2
+&&\Rightarrow&
+\Delta U &= -\lambda^2 U
+\end{align*}
+Die letzte Gleichung für die Funktion $U(x)$ hat die Form
+eines Eigenwertproblems mit dem Eigenwert $-\lambda^2$.
+
+\begin{definition}
+Eine Eigenfunktion eines Operators $L$ zum Eigenwert $\lambda$
+ist eine Funktion $U$ derart, dass $LU=\lambda U$.
+\end{definition}
+
+Die Separation ermöglich also, das ursprüngliche Problem aufzuspalten
+in ein Eigenwertproblem für eine nur ortsabhängige Funktion $U(x)$
+und eine Schwingungsgleichung für $T(t)$.
+Die Schwingungsfrequenz $c \lambda $ hängt direkt mit dem
+Eigenwert zusammen.
+Die Funktion $U(x)$ beschreibt die Form der Membran, die Amplitude
+in jedem Punkt, der Faktor $T(t)$ beschreibt die Schwingung.
+
+
diff --git a/buch/chapters/part1.tex b/buch/chapters/part1.tex
index b449aa2..51134ba 100644
--- a/buch/chapters/part1.tex
+++ b/buch/chapters/part1.tex
@@ -19,7 +19,7 @@
 \input{chapters/060-integral/chapter.tex}
 %\input{chapters/070-reihenprodukte/chapter.tex}
 \input{chapters/080-funktionentheorie/chapter.tex}
-%\input{chapters/090-funktional/chapter.tex}
+\input{chapters/090-pde/chapter.tex}
 % Gamma und Pi
 % Eulersche Beta-Funktion
 
diff --git a/buch/common/macros.tex b/buch/common/macros.tex
index 2c6eea2..7c82180 100644
--- a/buch/common/macros.tex
+++ b/buch/common/macros.tex
@@ -104,6 +104,7 @@
 \newtheorem{frage}[satz]{Frage}
 \newtheorem{problem}[satz]{Problem}
 \newtheorem{aufgabe}[satz]{Aufgabe}
+\newtheorem{prinzip}[satz]{Prinzip}
 \newtheorem*{problem*}{Problem}
 \newtheorem{forderung}{Forderung}[chapter]
 \newtheorem{konsequenz}[satz]{Konsequenz}