Crank-Nicolson-Methode

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

Crank–Nicolson beim 1D-Problem.

Das Crank-Nicolson-Verfahren basiert auf der Finite-Differenzen-Methode im Raum und 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)$

und sei $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 mit der Finite-Differenzen-Methode räumlich diskretisiert werden.

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

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

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

die Crank–Nicolson-Diskretisierung ist 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 Matrix-Inversion 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. 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.

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

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 in Folge 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., α, β und λ 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

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 Δx = Δ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

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.^[3] 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

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. ↑ The Mathematics of Financial Derivatives: A Student Introduction. Cambridge Univ. Press 1995, ISBN 0521497892

Weblinks

• Modul für parabolische partielle DGL (engl.)