v3.2 – Instationäre Diffusionsprobleme

 

3.2.1 Aufgabenstellung

Wir betrachten die Anfangs-Randwertaufgabe

\frac{{\partial u}}{{\partial t}}-\frac{{{\partial ^2}u}}{{\partial {x^2}}} = f\left( {x,t} \right)\quad \quad in\:\:\:\left( {0,1} \right) \times \left( {0,T} \right)

u\left( {x,0} \right) = {u_0}\left( x \right)

u\left( {0,t} \right) = u\left( {1,t} \right) = 0

3.2.2 Diskretisierung mit Differenzenverfahren

Wir definieren ein Ortsgitter der Schrittweite h = \frac{1}{N}

{x_j} = jh,\quad j = 0, \ldots ,N

und ein Zeitgitter mit der Schrittweite \tau = \frac{T}{M}

{t^k} = k\tau ,\quad k = 0, \ldots ,M

und approximieren die Lösung u durch die Gitterfunktion {u_{h,\tau }}, wobei wir wieder

u_j^k: = {u_{h,\tau }}\left( {{x_j},{t^k}} \right)

definieren. Die Zeitableitung wird standardmäßig durch einen Differenzenquotienten ersetzt:

\frac{{\partial u}}{{\partial t}} \approx \frac{1}{\tau }\left( {u_j^{k+1}-u_j^k} \right)

Die zweite partielle Ableitung im Ort wird durch eine zweite Differenz ersetzt:

\frac{{{\partial ^2}u}}{{\partial {x^2}}} \approx \frac{1}{{{h^2}}}\left( {{u_{j+1}}-2{u_j}+{u_{j-1}}} \right)

wobei sich die Frage stellt, auf welcher Zeitschicht die Gitterfunktion ausgewertet wird. Für maximale Flexibilität führt man ein Gewicht \sigma \in \left[ {0,1} \right] ein und approximiert

\frac{{{\partial ^2}u}}{{\partial {x^2}}} \approx \sigma \frac{1}{{{h^2}}}\left( {u_{j+1}^{k+1}-2u_j^{k+1}+u_{j-1}^{k+1}} \right)+\left( {1-\sigma } \right)\frac{1}{{{h^2}}}\left( {u_{j+1}^k-2u_j^k+u_{j-1}^k} \right)

Diese Approximation wird als \sigma-gewichtetes Schema bezeichnet. Für die Approximation der rechten Seite schreiben wir \tilde f_j^k und lassen zunächst noch offen, ob wir \tilde f_j^k = f\left( {{x_j},{t^k}} \right) oder \tilde f_j^k = f\left( {{x_j},{t^k}+\frac{\tau }{2}} \right) wählen. Die Diskretisierung lautet folglich:

\frac{1}{\tau }\left( {u_j^{k+1}-u_j^k} \right)-\sigma \frac{1}{{{h^2}}}\left( {u_{j+1}^{k+1}-2u_j^{k+1}+u_{j-1}^{k+1}} \right)

-\left( {1-\sigma } \right)\frac{1}{{{h^2}}}\left( {u_{j+1}^k-2u_j^k+u_{j-1}^k} \right) = \tilde f_j^k

mit j = 1, \ldots ,N-1 und k = 0, \ldots ,M-1. Die Anfangs- und Randbedingungen werden standardmäßig diskretisiert:

u_j^0 = {u_0}\left( {{x_j}} \right),\quad j = 0, \ldots ,N

u_0^k = u_N^k = 0,\quad k = 0, \ldots ,M

3.2.3 Differenzengleichung als Gleichungssystem

Wir nehmen nun an, dass alle Werte für {t^k} vorliegen und die Werte für {t^{k+1}} berechnet werden sollen. Die Differenzengleichung kann dann in Form des Gleichungssystems

-\sigma \frac{\tau }{{{h^2}}}u_{j+1}^{k+1}+\left( {1+2\sigma \frac{\tau }{{{h^2}}}} \right)u_j^{k+1}-\sigma \frac{\tau }{{{h^2}}}u_{j-1}^{k+1} = F_j^k

geschrieben werden, wobei

F_j^k: = \tau \tilde f_j^k+\left( {1-\sigma } \right)\frac{\tau }{{{h^2}}}u_{j+1}^k+\left( {1-\sigma } \right)\frac{\tau }{{{h^2}}}u_{j-1}^k+\left( {1-2\left( {1-\sigma } \right)\frac{\tau }{{{h^2}}}} \right)u_j^k

gesetzt wurde.

Für \sigma = 0 vereinfachen sich die Gleichungen zu u_j^{k+1} = F_j^k, wir erhalten ein explizites Verfahren.

Für \sigma > 0 ist ein trilineares Gleichungssystem zu lösen, man spricht von einem impliziten Verfahren. Zwei Werte von \sigma > 0 haben eine herausgehobene Bedeutung:

Das Verfahren mit \sigma = 1 bezeichnet man als rein implizites Verfahren. Es zeichnet sich durch besondere Stabilität aus.

Für \sigma = 0,5 ergibt sich das Crank-Nicolson-Verfahren, wobei

\tilde f_j^k = f\left( {{x_j},{t^k}+\frac{\tau }{2}} \right) oder \tilde f_j^k = \frac{1}{2}\left( {f\left( {{x_j},{t^k}} \right)+f\left( {{x_j},{t^{k+1}}} \right)} \right)

gewählt wird. Dies zeichnet sich durch eine höhere Konsistenz- bzw. Konvergenzordnung aus.

3.2.4 Konsistenzfehler

Für hinreichend glatte Lösungen ergibt sich ein Konsistenzfehler von O\left( {{\tau ^2}+{h^2}} \right) für das Crank-Nicolson-Verfahren und O\left( {\tau +{h^2}} \right) für \sigma \ne 0,5.

3.2.5 Stabilität

Stabilität kann in verschiedenen Normen untersucht werden. Es ergeben sich folgende Stabilitätsbedingungen:

  • Stabilität in der Maximumnorm: \left( {1-\sigma } \right)\tau \leq \frac{1}{2}{h^2}
  • Stabilität in der schwächeren Euklidischen Norm (Euklidische Norm multipliziert mit h, diskretisierte {L^2}-Norm): \left( {1-2\sigma } \right)\tau \leq \frac{1}{2}{h^2}

Das bedeutet für die drei wichtigsten Verfahren:

  • \sigma = 0: In beiden Normen die Stabilitätsbedingung \tau \leq \frac{1}{2}{h^2}
    Dies ist sehr klein, aber in Ordnung, da der Konsistenzfehler auch O\left( {\tau +{h^2}} \right)
  • \sigma = 1: unbedingte Stabilität in beiden Normen
  • \sigma = 0,5: unbedingte Stabilität in der Euklidischen Norm, aber Stabilität in der Maximumnorm nur, wenn \tau \leq {h^2}

3.2.6 Konvergenz

Aus Stabilität und Konsistenz folgt wieder Konvergenz.

Für das Crank-Nicolson-Verfahren vergleichen wir hier mit der Schrittweite h = 0,1 die beiden Varianten \tau = 0,01 und \tau = 0,07:

num-302-instationares-diffusionsproblem-maximumprinzip-konvergenz

Bei \tau = 0,07 ist das Verfahren nur in der Euklidischen Norm stabil. Es konvergiert zwar, aber die Maximumnorm ist verletzt. Das liegt daran, dass die Euklidische Norm Durchschnitte bildet, statt das Maximum zu suchen.

3.2.7 Teilschritt-Theta-Verfahren

Eine Anwendungsvariante des \sigma-gewichteten Schemas besteht darin, statt eines Zeitschritts \tau drei Teilschritte der Länge

{\tau _1} = \theta \tau ,\quad {\tau _2} = \left( {1-2\theta } \right)\tau ,\quad {\tau _3} = \theta \tau

mit \theta = 1-\frac{1}{2}\sqrt 2 \approx 0,3 und den Gewichten

{\sigma _1} = \alpha ,\quad {\sigma _2} = 1-\alpha ,\quad {\sigma _3} = \alpha

mit \alpha \in \left( {\frac{1}{2},1} \right] auszuführen. Man erhält das Teilschritt-Theta-Verfahren (fractional step theta method).

Dieses Verfahren ist besonders stabil (stark A-stabil) und ebenfalls von zweiter Ordnung, wobei die Fehlerkonstante nur wenig größer als beim Crank-Nicolson-Verfahren ist. Es ist besonders geeignet zur Behandlung von parabolischen Problemen mit nicht regulären Daten.

Wählt man \alpha = \frac{{1-2\theta }}{{1-\theta }} \approx 0,6, dann sind die Matrizen der in den 3 Teilschritten zu lösenden Gleichungssysteme gleich, was zum Beispiel bei der direkten Lösung ausgenutzt werden kann.

Ähnliche Artikel

Kommentar verfassen