v2.1 – Finite Differenzen für die Zwei-Punkt-Randwertaufgabe

 

Randwert- und Randeigenwertaufgaben (Finite Differenzen)

2.1.1 Wärmeleitgleichung

Wir suchen die Temperaturverteilung in einem Stab:

num-201-stab-temperatur-warmeleitung-differentialgleichung

Der Körper A hat die Temperatur {g_0} = const, der Körper B hat die Temperatur {g_1} = const.

Das mathematische Modell lautet:

\underbrace {-{{\left( {k\left( x \right){u^\prime }\left( x \right)} \right)}^\prime }+q\left( x \right)u\left( x \right)}_{ = Lu} = f\left( x \right)+q\left( x \right){u_0}\left( x \right)\quad \quad in\:\:\:\Omega = \left( {0,1} \right)

u\left( 0 \right) = {g_0},\quad \quad u\left( 1 \right) = {g_1}

Dabei bedeuten

u\left( x \right): Temperatur im Stab
{u_0}\left( x \right): Umgebungstemperatur
k\left( x \right) \geq {k_0} > 0: Wärmeleitfähigkeit
q\left( x \right) > 0: Wärmeaustauschkoeffizient
f\left( x \right): Dichte der im Stab vorhandenen Energiequellen
L: Differentialoperator 2. Ordnung

Wir haben eine gewöhnliche Differentialgleichung 2. Ordnung und je eine Bedingung an den Randpunkten des Rechengebietes \Omega, also eine Randwertaufgabe.

Handelt es sich nicht um einen Stab, sondern um einen allgemeinen zwei- oder dreidimensionalen Körper, dann erhält man

Lu: = -\nabla \cdot \left( {k\left( x \right)\nabla u\left( x \right)} \right)+q\left( x \right)u\left( x \right) = f\left( x \right)\quad \quad in\:\:\:\Omega

u\left( x \right) = g\left( x \right)\quad \quad auf\:\:\Gamma = \partial \Omega

Diese Aufgaben sind nur für einfache Gebiete (Kreis, Rechteck) geschlossen lösbar.

2.1.2 Maximumprinzip, Vergleichsprinzip

Für das Wärmeleitproblem setzen wir q\left( x \right) \geq 0 voraus. Dann gelten die folgenden beiden Prinzipien:

Maximumprinzip:
Wenn f\left( x \right) \leq 0\quad in\:\:\:\Omega und {g_0},{g_1} \leq 0, dann ist u\left( x \right) \leq 0\:\:\forall x \in \Omega. Das Maximum wird auf dem Rand eingenommen.

Vergleichsprinzip:
Gilt für zwei Funktionen v und w:

Lv \leq Lw\quad in\:\:\:\Omega ,\quad v\left( 0 \right) \leq w\left( 0 \right),\quad v\left( 1 \right) \leq w\left( 1 \right)

dann ist v\left( x \right) \leq w\left( x \right)\:\:\forall x \in \Omega. Diese Eigenschaft wird auch als inverse Monotonie von L bezeichnet.

2.1.3 Schießverfahren

Um die Methoden aus Kapitel 1 auf die Randwertaufgaben anwenden zu können, formen wir die Differentialgleichung 2. Ordnung in ein System von DGL erster Ordnung um:

-{u^{\prime \prime }}+qu = f,\quad \quad u\left( 0 \right) = {g_0},\quad u\left( 1 \right) = {g_1}

\quad \Rightarrow \quad {u^\prime } = v,\quad {v^\prime } = qu-f

Nun bestimmen wir Näherungslösungen \Omega = \left( {0,1} \right),\quad {\Omega _h} = \left\{ {{x_i} = ih,\:\:i = 1, \ldots ,N-1} \right\} mit den Anfangswerten

{u_\tau }\left( {0,\alpha } \right) = {g_0}

{u_\tau }\left( {0,\alpha } \right) = \alpha

wobei der unbekannte Parameter \alpha so bestimmt werden muss, dass {u_\tau }\left( {1,\alpha } \right) = {g_1} erfüllt wird. Man kann das Verfahren als Lösung der nichtlinearen Gleichung {u_\tau }\left( {1,\alpha } \right) = {g_1} interpretieren, wobei die Berechnung eines Funktionswertes {u_\tau }\left( {1,\alpha } \right) die numerische Lösung einer Anfangswertaufgabe beinhaltet.

2.1.4 Methode der finiten Differenzen

2.1.4.1 Definitionen

Einfache Verfahren zur Lösung von Randwertaufgaben sind Differenzenverfahren (Methode der finiten Differenzen, FDM). Dabei ersetzt man

  • das Gebiet \Omega durch ein Gitter {\Omega _h}, z.B. \Omega = \left( {0,1} \right),\quad {\Omega _h} = \left\{ {{x_i} = ih,\:\:i = 1, \ldots ,N-1} \right\}, mit der Schrittweite h = \frac{1}{N}, und den Rand \Gamma durch {\Gamma _h} hier im Beispiel

    \Gamma = \left\{ {0,1} \right\},\quad {\Gamma _h} = \left\{ {{x_0},{x_N}} \right\}:

    num-201-gitter-diskretisierung-eindimensional

    Analog zu \bar \Omega = \Omega \cap \Gamma definieren wir auch {\bar \Omega _h} = {\Omega _h} \cap {\Gamma _h}

  • die Funktionen der Veränderlichen x \in \Omega durch Gitterfunktionen, die nur auf {\Omega _h} definiert sind:

    u:\bar \Omega \to \mathbb{R}\quad \Rightarrow \quad {u_h}:{{\bar \Omega }_h} \to \mathbb{R}

    f:\Omega \to \mathbb{R}\quad \Rightarrow \quad {f_h} = {R_h}f:{\Omega _h} \to \mathbb{R}

    wobei {R_h} der Restriktionsoperator ist

  • Ableitungen (Differentialquotienten) durch Differenzenquotienten:

    {u^\prime }\left( x \right)\quad \Rightarrow \quad \left( {{D^+}u} \right)\left( x \right) = \frac{{u\left( {x+h} \right)-u\left( x \right)}}{h} (Vorwärtsdifferenz)

    oder: {u^\prime }\left( x \right)\quad \Rightarrow \quad \left( {{D^-}u} \right)\left( x \right) = \frac{{u\left( x \right)-u\left( {x-h} \right)}}{h} (Rückwärtsdifferenz)

    oder: {u^\prime }\left( x \right)\quad \Rightarrow \quad \left( {{D^0}u} \right)\left( x \right) = \frac{{u\left( {x+h} \right)-u\left( {x-h} \right)}}{{2h}} (zentrale Differenz)

    {u^{\prime \prime }}\left( x \right)\quad \Rightarrow \quad \left( {{D^+}{D^-}u} \right)\left( x \right) = \frac{{\frac{{u\left( {x+h} \right)-u\left( x \right)}}{h}-\frac{{u\left( x \right)-u\left( {x-h} \right)}}{h}}}{h}

    \quad \quad \quad \quad \quad \quad \quad \quad \quad = \frac{{u\left( {x+h} \right)-2u\left( x \right)+u\left( {x-h} \right)}}{{{h^2}}}

Man erhält die endlichdimensionale Ersatzaufgabe

{L_h}{u_h} = {f_h}\quad in\:\:\:{\Omega _h}

2.1.4.2 Beispiel

Wir betrachten die Randwertaufgabe

-{u^{\prime \prime }}+qu = f

u\left( 0 \right) = {g_0}

u\left( 1 \right) = {g_1}

Die Approximation lautet:

-{D^+}{D^-}{u_h}+{c_h}{u_h} = {f_h}

{u_h}\left( 0 \right) = {g_0}

{u_h}\left( 1 \right) = {g_1}

Für die Punkte {x_k},\:\:k = 1, \ldots ,N-1 bedeutet das:

-\frac{{{u_{k+1}}-2{u_k}+{u_{k-1}}}}{{{h^2}}}+{c_k}{u_k} = {f_k}

wobei {u_k}: = {u_h}\left( {{x_k}} \right),\:\:\:{c_k}: = {c_h}\left( {{x_k}} \right) = c\left( {{x_k}} \right),\:\:\:{f_k}: = {f_h}\left( {{x_k}} \right) = f\left( {{x_k}} \right) gesetzt wurde.

Die Approximation kann auch als lineares Gleichungssystem geschrieben werden:

\left( {\begin{array}{*{20}{c}}  1 & {} & {} & {} & {} & {} \\{-\frac{1}{{{h^2}}}} & {{c_1}+\frac{2}{{{h^2}}}} & {-\frac{1}{{{h^2}}}} & {} & {} & {} \\{} & {-\frac{1}{{{h^2}}}} & {{c_2}+\frac{2}{{{h^2}}}} & {-\frac{1}{{{h^2}}}} & {} & {} \\{} & {} & \ddots & \ddots & \ddots & {} \\{} & {} & {} & {-\frac{1}{{{h^2}}}} & {{c_{N-1}}+\frac{2}{{{h^2}}}} & {-\frac{1}{{{h^2}}}} \\{} & {} & {} & {} & {} & 1 \\   \end{array} } \right)\left( {\begin{array}{*{20}{c}}{{u_0}} \\{{u_1}} \\{{u_2}} \\  \vdots \\{{u_{N-1}}} \\{{u_N}} \\   \end{array} } \right) = \left( {\begin{array}{*{20}{c}}{{g_0}} \\{{f_1}} \\{{f_2}} \\  \vdots \\{{f_{N-1}}} \\{{g_1}} \\   \end{array} } \right)

Die Systemmatrix ist großdimensioniert, schwach besetzt (nur O\left( N \right) Nicht-Nullelemente) und schlecht konditioniert (Konditionszahl O\left( {{N^2}} \right)). Oft eliminiert man die erste und letzte Gleichung.

2.1.5 Matrixnorm und Konditionszahl

Zeilensummennorm: {\left\| A \right\|_\infty } = \max \limits_{i = 1, \ldots ,n} \sum\limits_{j = 1}^m {\left| {{a_{ij}}} \right|}

Spektralnorm: {\left\| A \right\|_2} = \sqrt {{\lambda _{\max} }\left( {{A^T}A} \right)}, bei symmetrischer Matrix: {\left\| A \right\|_2} = {\lambda _{\max} }\left( A \right)

Frobeniusnorm: {\left\| A \right\|_F} = \sqrt {\sum\limits_{i,j = 1}^n {{{\left| {{a_{ij}}} \right|}^2}} }

Bedingung für verträgliche Norm: \left\| {Ax} \right\| \leq \left\| A \right\|\left\| x \right\|

Dies bekommt man mit induzierter Matrixnorm: \left\| A \right\|: = \sup \limits_{x \ne 0} \frac{{\left\| {Ax} \right\|}}{{\left\| x \right\|}} = \sup \limits_{\left\| x \right\| = 1} \left\| {Ax} \right\|

Konditionszahl: \kappa = \left\| A \right\|\left\| {{A^{-1}}} \right\|

Also unter Verwendung der Spektralnorm:

{\kappa _2} = {\left\| A \right\|_2}{\left\| {{A^{-1}}} \right\|_2} =  {\lambda _{\max} }\left( A \right){\lambda _{\max} }\left( {{A^{-1}}} \right) = \frac{{{\lambda _{\max} }\left( A \right)}}{{{\lambda _{\min} }\left( A \right)}}

2.1.6 Diskretes Maximumprinzip, diskretes Vergleichsprinzip

Diskretes Maximumprinzip:
Wenn {f_i} \leq 0,\:\:i = 1, \ldots ,N-1,\:\:\:{g_0} \leq 0,\:\:{g_1} \leq 0\quad \Rightarrow \quad {u_i} \leq 0\:\:\forall i = 0, \ldots ,N

Diskretes Vergleichsprinzip:
Wenn {L_h}{v_h} \leq {L_h}{w_h}\:\:\:in\:\:\Omega ,\quad {v_0} \leq {w_0},\quad {v_N} \leq {w_N}\quad \Rightarrow \quad {v_h} \leq {w_h}\:\:\:in\:\:{\Omega _h}

Beweis diskretes Maximumprinzip:
Sei j derjenige Index, für den {u_j} =  \max \limits_i {u_i}. Wenn j = 0 oder j = N ist, dann ist alles gezeigt. Sei also 0 < j < N. Dann gilt {u_{j+1}}-{u_j} \leq 0,\quad {u_{j-1}}-{u_j} \leq 0. Daraus folgt:

\underbrace {{c_j}}_{ > 0}{u_j} = {f_j}+\frac{{{u_{j+1}}-2{u_j}+{u_{j-1}}}}{{{h^2}}} = \underbrace {{f_j}}_{ \leq 0}+\underbrace {\frac{{{u_{j+1}}-{u_j}}}{{{h^2}}}}_{ \leq 0}+\underbrace {\frac{{{u_{j-1}}-{u_j}}}{{{h^2}}}}_{ \leq 0} \leq 0\quad \Rightarrow \quad {u_j} \leq 0

Das diskrete Vergleichsprinzip folgt aus dem Maximumprinzip mit {u_h} = {v_h}-{w_h}.

Aus dem diskreten Vergleichsprinzip folgt, dass die Aufgabe

{L_h}{u_h} = 0\quad in\:\:\:{\Omega _h},\quad {u_h}\left( 0 \right) = 0,\quad {u_h}\left( 1 \right) = 0

nur die Null-Lösung besitzt. Daher besitzt die Diskretisierung

{L_h}{u_h} = f\quad in\:\:\:{\Omega _h},\quad {u_h}\left( 0 \right) = {g_0},\quad {u_h}\left( 1 \right) = {g_1}

eine eindeutige Lösung.

2.1.7 Diskrete Konvergenz, Konsistenz und Stabilität

2.1.7.1. Definitionen

Für Gitterfunktionen sei die Maximumnorm durch

{\left\| {{v_h}} \right\|_{{\Omega _h}}} = \max \limits_{{x_i} \in {\Omega _h}} \left| {{v_h}\left( {{x_i}} \right)} \right|,\quad \quad {\left\| {{v_h}} \right\|_{{{\bar \Omega }_h}}} = \max \limits_{{x_i} \in {{\bar \Omega }_h}} \left| {{v_h}\left( {{x_i}} \right)} \right|

definiert. Ein Differenzenverfahren heißt diskret konvergent mit der Ordnung p, wenn

{\left\| {{R_h}u-{u_h}} \right\|_{{{\bar \Omega }_h}}} = O\left( {{h^p}} \right)

Ein Differenzenverfahren heißt konsistent mit der Ordnung p, wenn

{\left\| {{L_h}{R_h}u-{R_h}Lu} \right\|_{{{\bar \Omega }_h}}} = O\left( {{h^p}} \right)

Ein Differenzenverfahren heißt stabil bezüglich der rechten Seite, wenn eine Konstante {C_S} existiert, so dass für alle Gitterfunktionen {v_h} mit {v_h} = 0 auf {\Gamma _h} die Abschätzung

{\left\| {{v_h}} \right\|_{{{\bar \Omega }_h}}} \leq {C_S}{\left\| {{L_h}{v_h}} \right\|_{{\Omega _h}}}

gilt.

Ist ein Differenzenverfahren stabil und konsistent mit der Ordnung p, dann ist es auch konvergent mit der Ordnung p.

Beweis:

Es gilt

{\left\| {{R_h}u-{u_h}} \right\|_{{{\bar \Omega }_h}}} \leq {C_S}{\left\| {{L_h}\left( {{R_h}u-{u_h}} \right)} \right\|_{{\Omega _h}}} = {C_S}{\left\| {{L_h}{R_h}u-{L_h}{u_h}} \right\|_{{\Omega _h}}}

Wegen {L_h}{u_h} = {f_h} = {R_h}f = {R_h}Lu folgt

{\left\| {\underbrace {{R_h}u-{u_h}}_{Diskretisierungsfehler}} \right\|_{{{\bar \Omega }_h}}} \leq {C_S}{\left\| {\underbrace {{L_h}{R_h}u-{R_h}Lu}_{Konsistenzfehler}} \right\|_{{\Omega _h}}}

Je kleiner {C_S}, desto besser.

2.1.7.2 Untersuchung von Konsistenz

Konsistenz untersucht man über Taylorreihenentwicklung. Im Punkt {x_i} gilt für unsere Beispielaufgabe

\left( {{L_h}{R_h}u-{R_h}Lu} \right)\left( {{x_i}} \right)

= -\frac{{u\left( {{x_i}+h} \right)-2u\left( {{x_i}} \right)+u\left( {{x_i}-h} \right)}}{{{h^2}}}+c\left( {{x_i}} \right)u\left( {{x_i}} \right)-\left( {-{u^{\prime \prime }}\left( {{x_i}} \right)+c\left( {{x_i}} \right)u\left( {{x_i}} \right)} \right)

= {u^{\prime \prime }}\left( {{x_i}} \right)-\frac{{u\left( {{x_i}+h} \right)-2u\left( {{x_i}} \right)+u\left( {{x_i}-h} \right)}}{{{h^2}}}

Es gilt nun für ein mindestens viermal differenzierbares u:

u\left( {{x_i} \pm h} \right) = u\left( {{x_i}} \right) \pm h{u^\prime }\left( {{x_i}} \right)+\frac{{{h^2}}}{2}{u^{\prime \prime }}\left( {{x_i}} \right) \pm \frac{{{h^3}}}{6}{u^{\left( 3 \right)}}\left( {{x_i}} \right)+\frac{{{h^4}}}{{4!}}{u^{\left( 4 \right)}}\left( {{\xi _ \pm }} \right)

{\xi _-} \in \left( {{x_i}-h,{x_i}} \right),\quad {\xi _+} \in \left( {{x_i},{x_i}+h} \right)

Eingesetzt:

u\left( {{x_i}+h} \right)+u\left( {{x_i}-h} \right) = 2u\left( {{x_i}} \right)+{h^2}{u^{\prime \prime }}\left( {{x_i}} \right)+O\left( {{h^4}} \right)

\quad \Rightarrow \quad \left( {{L_h}{R_h}u-{R_h}Lu} \right)\left( {{x_i}} \right) = {u^{\prime \prime }}\left( {{x_i}} \right)-\frac{{u\left( {{x_i}+h} \right)-2u\left( {{x_i}} \right)+u\left( {{x_i}-h} \right)}}{{{h^2}}}

\quad \Rightarrow \quad \left( {{L_h}{R_h}u-{R_h}Lu} \right)\left( {{x_i}} \right) = {u^{\prime \prime }}\left( {{x_i}} \right)-\frac{{2u\left( {{x_i}} \right)+{h^2}{u^{\prime \prime }}\left( {{x_i}} \right)+O\left( {{h^4}} \right)-2u\left( {{x_i}} \right)}}{{{h^2}}}

\quad \Rightarrow \quad \left( {{L_h}{R_h}u-{R_h}Lu} \right)\left( {{x_i}} \right) = O\left( {{h^2}} \right)

\quad \Rightarrow \quad {\left\| {{L_h}{R_h}u-{R_h}Lu} \right\|_{{{\bar \Omega }_h}}} = O\left( {{h^2}} \right)

Das Verfahren ist demnach konsistent mit Ordnung 2.

2.1.7.3 Untersuchung von Stabilität

Wir zeigen die Stabilität in unserer Beispielaufgabe, wobei wir zur Vereinfachung c\left( x \right) \geq {c_0} > 0 voraussetzen. Sei {v_h} eine beliebige Gitterfunktion mit {v_h} = 0 auf {\Gamma _h}. Ist {L_h}{v_h} \equiv 0, dann ist {v_h} \equiv 0 und wir sind fertig. Sei nun {L_h}{v_h} \ne 0. Wir betrachten die Hilfsfunktion

{w_h} = \frac{{{{\left\| {{L_h}{v_h}} \right\|}_{{\Omega _h}}}}}{{{c_0}}} > 0

Da {w_h} konstant ist, gilt für alle {x_i} \in {\Omega _h}:

\left( {{L_h}{w_h}} \right)\left( {{x_i}} \right) = c\left( {{x_i}} \right){w_i} \geq {c_0}{w_i} = {\left\| {{L_h}{v_h}} \right\|_{{\Omega _h}}} \geq \pm {L_h}{v_h}\left( {{x_i}} \right)

also

-{L_h}{w_h}\left( {{x_i}} \right) \leq {L_h}{v_h}\left( {{x_i}} \right) \leq {L_h}{w_h}\left( {{x_i}} \right)

Da andererseits -{w_h} < 0 = {v_h} < {w_h} auf {\Gamma _h} gilt, folgt mit dem diskreten Vergleichsprinzip -{w_h} \leq {v_h} \leq {w_h} in {\bar \Omega _h} und damit

{\left\| {{v_h}} \right\|_{{{\bar \Omega }_h}}} \leq {\left\| {{w_h}} \right\|_{{{\bar \Omega }_h}}} = \frac{{{{\left\| {{L_h}{v_h}} \right\|}_{{\Omega _h}}}}}{{{c_0}}}

also Stabilität mit der Konstanten {C_S} = c_0^{-1}.

Insgesamt haben wir also Konvergenz mit Ordnung 2.