Crank-Nicolson-Verfahren

aus Wikipedia, der freien Enzyklopädie
Wechseln zu: Navigation, Suche

Das Crank-Nicolson-Verfahren ist in der numerischen Mathematik eine Finite-Differenzen-Methode zur Lösung der Wärmeleitungsgleichung und ähnlicher partieller Differentialgleichungen.[1] Es ist ein implizites Verfahren 2. Ordnung und numerisch stabil. Das Verfahren wurde Mitte des 20. Jahrhunderts von John Crank und Phyllis Nicolson entwickelt.[2]

Für die Wärmeleitungsgleichung und viele andere Gleichungen kann gezeigt werden, dass das Crank-Nicolson-Verfahren ohne Bedingungen numerisch stabil ist. Trotzdem können die approximierten Lösungen störende Schwingungen enthalten, wenn der Quotient aus Zeitdifferenz und Abstandsquadrat groß ist (typischerweise größer als 1/2). In diesem Fall wird häufig das weniger genaue Euler-Rückwärtsverfahren genutzt, welches numerisch stabil und unempfindlich gegenüber Störungen ist.

Das Verfahren[Bearbeiten]

Crank–Nicolson beim 1D-Problem.

Das Crank-Nicolson-Verfahren basiert auf der Trapezregel in der Zeit. Dies ist gleichbedeutend mit dem Mittelwert des Euler-Vorwärtsverfahren und des Euler-Rückwärtsverfahren.

Ist im eindimensionalen Fall die partielle Differentialgleichung

\frac{\partial u}{\partial t} = F\left(u, x, t, \frac{\partial u}{\partial x}, \frac{\partial^2 u}{\partial x^2}\right)

gegeben und bezeichne u(i \Delta x, n \Delta t) = u_{i}^{n}\,, dann ist das Crank-Nicolson-Verfahren der Mittelwert des Euler-Vorwärtsverfahren beim Wert n und des Euler-Rückwärtsverfahren bei n + 1:

\frac{u_{i}^{n + 1} - u_{i}^{n}}{\Delta t} = 
F_{i}^{n}\left(u, x, t, \frac{\partial u}{\partial x}, \frac{\partial^2 u}{\partial x^2}\right) \qquad \mbox{(Euler-Vorwärts)}
\frac{u_{i}^{n + 1} - u_{i}^{n}}{\Delta t} = 
F_{i}^{n + 1}\left(u, x, t, \frac{\partial u}{\partial x}, \frac{\partial^2 u}{\partial x^2}\right) \qquad \mbox{(Euler-Rückwärts)}
\frac{u_{i}^{n + 1} - u_{i}^{n}}{\Delta t} = 
\frac{1}{2}\left(
F_{i}^{n + 1}\left(u, x, t, \frac{\partial u}{\partial x}, \frac{\partial^2 u}{\partial x^2}\right) + 
F_{i}^{n}\left(u, x, t, \frac{\partial u}{\partial x}, \frac{\partial^2 u}{\partial x^2}\right)
\right) \qquad \mbox{(Crank-Nicolson)}

Die Funktion F muss räumlich diskretisiert werden, z.B. mit der Finite-Differenzen-Methode, dem Finite-Volumen-Verfahren oder dem Finite-Elemente-Verfahren.

Es ist zu beachten, dass es sich um eine implizite Methode handelt. Um den zeitlich „nächsten“ Wert von u zu erhalten, muss ein algebraisches Gleichungssystem gelöst werden. Ist die partielle Differentialgleichung nicht linear, wird die Diskretisierung ebenfalls nicht linear sein, so dass ein nicht lineares algebraisches Gleichungssystem gelöst werden muss, obwohl Linearisierungen möglich sind. Bei vielen Problemen, insbesondere bei linearer Diffusion, ist das algebraische Problem tridiagonal, das am besten mit dem schnellen Thomas-Algorithmus gelöst wird, um eine aufwändige Inversion der Matrix zu vermeiden.

Beispiel: Eindimensionale Diffusion[Bearbeiten]

Das Crank-Nicolson-Verfahren wird häufig zur Lösung von Diffusionsproblemen verwendet. Ein Beispiel für lineare Diffusion ist

\frac{\partial u}{\partial t} = a \frac{\partial^2 u}{\partial x^2}.

Unter Verwendung der Finite-Differenzen-Methode für die räumliche Diskretisierung der rechten Seite, ist die Crank–Nicolson-Diskretisierung dann:

\frac{u_{i}^{n + 1} - u_{i}^{n}}{\Delta t} = \frac{a}{2 (\Delta x)^2}\left(
(u_{i + 1}^{n + 1} - 2 u_{i}^{n + 1} + u_{i - 1}^{n + 1}) + 
(u_{i + 1}^{n} - 2 u_{i}^{n} + u_{i - 1}^{n})
\right)

oder, sei r = \frac{a \Delta t}{2 (\Delta x)^2}:

-r u_{i + 1}^{n + 1} + (1 + 2 r)u_{i}^{n + 1} - r u_{i - 1}^{n + 1} = r u_{i + 1}^{n} + (1 - 2 r)u_{i}^{n} + r u_{i - 1}^{n}\,

dies ist ein tridiagonales Problem, so dass u_{i}^{n + 1}\, mit dem Thomas-Algorithmus gelöst werden sollte, um eine aufwändige Inversion der entstehenden Tridiagonal-Toeplitz-Matrix zu vermeiden.

Eine quasilineare Gleichung, wie (dies ist ein vereinfachtes Beispiel und nicht allgemein gültig)

\frac{\partial u}{\partial t} = a(u) \frac{\partial^2 u}{\partial x^2}

würde zu einem nicht linearen System von algebraischen Gleichungen führen, die nicht so einfach zu lösen sind wie oben. Es ist jedoch in manchen Fällen möglich, das Problem zu linearisieren, indem man den alten Wert von a nutzt; dieses Vorgehen wird auch als "Lagging the Coefficients" bezeichnet. Dieser ist a_{i}^{n}(u)\, statt a_{i}^{n + 1}(u)\,. In anderen Fällen kann es möglich sein, a_{i}^{n + 1}(u)\, mit Hilfe einer expliziten Methode und Beibehaltung der Stabilität zu schätzen.

Alternativ ermöglicht iteratives Aktualisieren der Koeffizienten, die Ordnung des Verfahren beizubehalten. Hierbei wird "Lagging the Coefficients" in einem ersten Schritt angewendet, um das Gleichungssystem für den nächsten Zeitschritt zu lösen. Mit den nun bekannten neuen Temperaturen werden die Koeffizienten ausgewertet und die Lösung wird erneut berechnet. Dieses Vorgehen wird fortgesetzt bis die Differenz der Lösungen zweier aufeinanderfolgender Iterationen gering ist. [3]

Beispiel: Eindimensionale Diffusion mit Advektion für stationäre Strömungen[Bearbeiten]

Diese Lösung wird gewöhnlich für Zwecke genutzt, wenn eine Kontamination in Fließgewässern mit stationärer Strömung vorliegt, aber nur Informationen in einer Dimension gegeben sind. Oft kann das Problem zu einem eindimensionalen Problem vereinfacht werden und dennoch nützliche Informationen ergeben.

Im Folgenden wird ein aufgelöster Fremdkörper in Wasser modelliert. Dieses Problem ist aus drei Teilen zusammengesetzt: die bekannte Diffusionsgleichung (gewählt als Konstante D_x), eine advektive Komponente (das System entwickelt sich infolge eines Vektorfeldes im Raum), welche als Konstante U_x gewählt wird, und eine laterale Überlagerung zwischen longitudinalen Kanälen (k)

\frac{\partial C}{\partial t} = D_x \frac{\partial^2 C}{\partial x^2} - U_x \frac{\partial C}{\partial x}- k (C-C_N)-k(C-C_M) \qquad\qquad (*)

wobei C die Konzentration des Fremdkörpers im Wasser ist und die Indizes N und M dem vorherigen und dem nächsten Kanal entsprechen.

Das Crank-Nicolson-Verfahren (wobei i für den Ort und j für die Zeit steht) transformiert jede Komponente der partiellen Differentialgleichung folgendermaßen:

1.\quad\frac{\partial C}{\partial t} = \frac{C_{i}^{j + 1} - C_{i}^{j}}{\Delta t}
2.\quad\frac{\partial^2 C}{\partial x^2}= \frac{1}{2 (\Delta x)^2}\left(
(C_{i + 1}^{j + 1} - 2 C_{i}^{j + 1} + C_{i - 1}^{j + 1}) + 
(C_{i + 1}^{j} - 2 C_{i}^{j} + C_{i - 1}^{j})
\right)
3.\quad\frac{\partial C}{\partial x} = \frac{1}{2}\left(
\frac{(C_{i + 1}^{j + 1} - C_{i - 1}^{j + 1})}{2 (\Delta x)} + 
 \frac{(C_{i + 1}^{j} - C_{i - 1}^{j})}{2 (\Delta x)}
\right)
4.\quad C= \frac{1}{2} (C_{i}^{j+1} + C_{i}^{j})
5.\quad C_N= \frac{1}{2} (C_{Ni}^{j+1} + C_{Ni}^{j})
6.\quad C_M= \frac{1}{2} (C_{Mi}^{j+1} + C_{Mi}^{j})

Nun werden folgende Konstanten definiert, um die Rechnung zu vereinfachen:

 \lambda= \frac{D_x\Delta t}{2 \Delta x^2}
 \alpha= \frac{U_x\Delta t}{4 \Delta x}
 \beta= \frac{k\Delta t}{2}

Nun werden 1. bis 6., \alpha, \beta und \lambda in (*) eingesetzt. Im Anschluss werden die „in der Zukunft liegenden“ Terme (j+1) auf die linke Seite und die „aktuellen“ (j) auf die rechte Seite gebracht. Es ergibt sich:

 -\beta C_{Ni}^{j+1}-(\lambda+\alpha)C_{i-1}^{j+1} +(1+2\lambda+2\beta)C_{i}^{j+1}-(\lambda-\alpha)C_{i+1}^{j+1}-\beta C_{Mi}^{j+1} = \beta C_{Ni}^{j}+(\lambda+\alpha)C_{i-1}^{j} +(1-2\lambda-2\beta)C_{i}^{j}+(\lambda-\alpha)C_{i+1}^{j}+\beta C_{Mi}^{j}

Um den ersten Kanal zu modellieren, wird festgestellt, dass es nur im Kontakt mit dem folgenden Kanal (M) sein kann, so dass sich der Ausdruck vereinfacht zu:

 -(\lambda+\alpha)C_{i-1}^{j+1} +(1+2\lambda+\beta)C_{i}^{j+1}-(\lambda-\alpha)C_{i+1}^{j+1}-\beta C_{Mi}^{j+1} = +(\lambda+\alpha)C_{i-1}^{j} +(1-2\lambda-\beta)C_{i}^{j}+(\lambda-\alpha)C_{i+1}^{j}+\beta C_{Mi}^{j}

Auf die gleiche Art und Weise wird der letzte Kanal modelliert. Es wird festgestellt, dass er nur im Kontakt mit dem vorherigen Kanal (N) sein kann, so dass sich der Ausdruck vereinfacht zu:

 -\beta C_{Ni}^{j+1}-(\lambda+\alpha)C_{i-1}^{j+1} +(1+2\lambda+\beta)C_{i}^{j+1}-(\lambda-\alpha)C_{i+1}^{j+1}= \beta C_{Ni}^{j}+(\lambda+\alpha)C_{i-1}^{j} +(1-2\lambda-\beta)C_{i}^{j}+(\lambda-\alpha)C_{i+1}^{j}

Um dieses lineare Gleichungssystem zu lösen, müssen Grenzbedingungen am Anfang der Kanäle gegeben sein:

 C_0^{j}: Anfangswert für den aktuellen Kanal
 C_{0}^{j+1}: Anfangswert für den Kanal des nächsten Zeitintervalls
 C_{N0}^{j}: Anfangswert für den vorherigen Kanal des aktuellen Zeitintervalls
 C_{M0}^{j}: Anfangswert für den nachfolgenden Kanal des aktuellen Zeitintervalls

Für die letzte Zelle der Kanäle (z) wird die günstigste Bedingung adiabatisch, also

\frac{\partial C}{\partial x}_{x=z} = \frac{(C_{i + 1} - C_{i - 1})}{2 \Delta x}  = 0

Diese Bedingung ist erfüllt, wenn (und nur dann (abgesehen vom Wert 0)):

 C_{i + 1}^{j+1} = C_{i - 1}^{j+1}

Es soll nun das Problem (in Matrizenform) für den Fall mit drei Kanälen und fünf Knoten (inklusive der Anfangswertbedingungen) gelöst werden. Als lineares Problem ausgedrückt:

 \begin{bmatrix}AA\end{bmatrix}\begin{bmatrix}C^{j+1}\end{bmatrix}=[BB][C^{j}]+[d]

wobei

\mathbf{C^{j+1}} = \begin{bmatrix}
C_{11}^{j+1}\\ C_{12}^{j+1} \\ C_{13}^{j+1} \\ C_{14}^{j+1} 
\\ C_{21}^{j+1}\\ C_{22}^{j+1} \\ C_{23}^{j+1} \\ C_{24}^{j+1} 
\\ C_{31}^{j+1}\\ C_{32}^{j+1} \\ C_{33}^{j+1} \\ C_{34}^{j+1} 
\end{bmatrix} \text{ und } \mathbf{C^{j}} = \begin{bmatrix}
C_{11}^{j}\\ C_{12}^{j} \\ C_{13}^{j} \\ C_{14}^{j} 
\\ C_{21}^{j}\\ C_{22}^{j} \\ C_{23}^{j} \\ C_{24}^{j} 
\\ C_{31}^{j}\\ C_{32}^{j} \\ C_{33}^{j} \\ C_{34}^{j} 
\end{bmatrix}

Es wird festgestellt, dass AA und BB Felder mit vier verschiedenen Feldkomponenten sein sollten. (Zur Erinnerung: es wurden für dieses Beispiel nur drei Kanäle berücksichtigt, aber es deckt dennoch den Hauptteil des oben Genannten ab)

\mathbf{AA} = \begin{bmatrix}
AA1 & AA3 & 0\\
AA3 & AA2 & AA3\\
0 & AA3 & AA1\end{bmatrix} \text{ und } \mathbf{BB} = \begin{bmatrix}
BB1 & -AA3 & 0\\
-AA3 & BB2 & -AA3\\
0 & -AA3 & BB1\end{bmatrix}  

wobei die oben erwähnten Elemente den nächsten Feldern und einer zusätzlichen 4x4-Nullmatrix entsprechen. Es ist zu beachten, dass AA und BB jeweils die Größe 12x12 haben:

\mathbf{AA1} = \begin{bmatrix}
(1+2\lambda+\beta) & -(\lambda-\alpha) & 0 & 0 \\
-(\lambda+\alpha) & (1+2\lambda+\beta) & -(\lambda-\alpha) & 0 \\
0 & -(\lambda+\alpha) & (1+2\lambda+\beta) & -(\lambda-\alpha)\\
0 & 0 & -2\lambda & (1+2\lambda+\beta)\end{bmatrix},
\mathbf{AA2} = \begin{bmatrix}
(1+2\lambda+2\beta) & -(\lambda-\alpha) & 0 & 0 \\
-(\lambda+\alpha) & (1+2\lambda+2\beta) & -(\lambda-\alpha) & 0 \\
0 & -(\lambda+\alpha) & (1+2\lambda+2\beta) & -(\lambda-\alpha)\\
0 & 0 & -2\lambda & (1+2\lambda+2\beta) \end{bmatrix},
\mathbf{AA3} = \begin{bmatrix}
-\beta & 0 & 0 & 0 \\
0 & -\beta & 0 & 0 \\
0 & 0 & -\beta & 0 \\
0 & 0 & 0 & -\beta \end{bmatrix},
\mathbf{BB1} = \begin{bmatrix}
(1-2\lambda-\beta) & (\lambda-\alpha) & 0 & 0 \\
(\lambda+\alpha) & (1-2\lambda-\beta) & (\lambda-\alpha) & 0 \\
0 & (\lambda+\alpha) & (1-2\lambda-\beta) & (\lambda-\alpha)\\
0 & 0 & 2\lambda & (1-2\lambda-\beta)\end{bmatrix}\text{ und}
\mathbf{BB2} = \begin{bmatrix}
(1-2\lambda-2\beta) & (\lambda-\alpha) & 0 & 0 \\
(\lambda+\alpha) & (1-2\lambda-2\beta) & (\lambda-\alpha) & 0 \\
0 & (\lambda+\alpha) & (1-2\lambda-2\beta) & (\lambda-\alpha)\\
0 & 0 & 2\lambda & (1-2\lambda-2\beta) \end{bmatrix}.

Der d-Vektor wird benutzt, um die Grenzbedingungen einzuhalten. In diesem Beispiel ist er ein 12x1-Vektor:

\mathbf{d} = \begin{bmatrix}
(\lambda+\alpha)(C_{10}^{j+1}+C_{10}^{j}) \\ 0 \\ 0 \\ 0 \\ (\lambda+\alpha)(C_{20}^{j+1}+C_{20}^{j}) \\ 0 \\ 0 \\ 0 \\ (\lambda+\alpha)(C_{30}^{j+1}+C_{30}^{j}) \\
0\\
0\\
0\end{bmatrix}

Um die Konzentration zu jedem Zeitpunkt zu finden, muss die folgende Gleichung iterativ gelöst werden:

 \left[ C^{j+1}\right]=\left[AA^{-1}\right]\left(\left[BB\right]\left[C^j\right]+\left[d\right]\right)

Beispiel: Zweidimensionale Diffusion[Bearbeiten]

Werden zwei Dimensionen betrachtet, ist die Herleitung analog und die Ergebnisse liefern eher ein System von band-diagonalen anstelle von tridiagonalen Gleichungen. Die zweidimensionale Wärmeleitungsgleichung

\frac{\partial u}{\partial t} = a \left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right)

kann gelöst werden mit der Crank-Nicolson-Diskretisierung

\begin{align}u_{i,j}^{n+1} &= u_{i,j}^n + \frac{1}{2} \frac{a \Delta t}{(\Delta x)^2} \big[(u_{i+1,j}^{n+1} + u_{i-1,j}^{n+1} + u_{i,j+1}^{n+1} + u_{i,j-1}^{n+1} - 4u_{i,j}^{n+1}) \\ & \qquad {} + (u_{i+1,j}^{n} + u_{i-1,j}^{n} + u_{i,j+1}^{n} + u_{i,j-1}^{n} - 4u_{i,j}^{n})\big]\end{align}

angenommen, dass ein rechtwinkliges Koordinatennetz genutzt wird, so dass \Delta x = \Delta y. Diese Gleichung kann ein wenig mit Umformungen des Terms und Benutzen der CFL-Zahl vereinfacht werden.

\mu = \frac{a \Delta t}{(\Delta x)^2}

Für die Stabilität des numerische Schemas von Crank-Nicolson ist eine niedrige CFL-Zahl nicht nötig, jedoch braucht man sie für die numerische Genauigkeit. Das Schema kann nun dargestellt werden als:

\begin{align}&(1 + 2\mu)u_{i,j}^{n+1} - \frac{\mu}{2}\left(u_{i+1,j}^{n+1} + u_{i-1,j}^{n+1} + u_{i,j+1}^{n+1} + u_{i,j-1}^{n+1}\right) \\ & \quad = (1 - 2\mu)u_{i,j}^{n} + \frac{\mu}{2}\left(u_{i+1,j}^{n} + u_{i-1,j}^{n} + u_{i,j+1}^{n} + u_{i,j-1}^{n}\right)\end{align}

Anwendung in der Finanzmathematik[Bearbeiten]

Weil auch eine Anzahl anderer Phänomene mit Hilfe der Wärmeleitungsgleichung (in der Finanzmathematik häufig Diffusionsgleichung genannt) modelliert werden können, wird das Crank-Nicolson-Verfahren auch in diesen Bereichen genutzt.[4] Insbesondere die Differentialgleichungen des Black-Scholes-Modells können in die Diffusionsgleichung umgewandelt und mit dem Crank-Nicolson-Verfahren gelöst werden. Die Wichtigkeit dessen kommt von der Ausbreitung des Optionspreismodells, das nicht mit einer in sich geschlossenen analytischen Lösung dargestellt werden kann, es können sich immer noch numerische Lösungen bieten. Allerdings ist das Crank-Nicolson-Verfahren bei ungleichmäßigen finanziellen Bedingungen (die bei den meisten Finanzinstrumenten vorliegen) nicht befriedigend, da numerische Schwingungen nicht gedämpft werden. Daher sind spezielle Dämpfungsinitialisierungen notwendig.

Einzelnachweise[Bearbeiten]

  1. Tuncer Cebeci: Convective Heat Transfer. Springer, 2002, ISBN 0966846141.
  2. J. Crank, P. Nicolson: A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. Springer, 1947.
  3. Pletcher, Richard H., John C. Tannehill, and Dale Anderson. Computational fluid mechanics and heat transfer. CRC Press, 2012.
  4. The Mathematics of Financial Derivatives: A Student Introduction. Cambridge Univ. Press, 1995, ISBN 0521497892.

Weblinks[Bearbeiten]