Crank-Nicolson-Verfahren

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

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

Differenzenstern (engl. stencil) für das Crank-Nicolson-Verfahren beim 1D-Problem. j ist der Ort und n ist die Zeit.

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

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ärtsverfahrens beim Wert $n$ und des Euler-Rückwärtsverfahrens 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 {\text{(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 {\text{(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 {\text{(Crank-Nicolson)}}$

Dabei stehen die Indizes $i$ und $n$ für den Ort bzw. die Zeit. 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 und wird am besten mit dem schnellen Thomas-Algorithmus gelöst, 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. 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}-2u_{i}^{n+1}+u_{i-1}^{n+1})+(u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n})\right)$

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

$-ru_{i+1}^{n+1}+(1+2r)u_{i}^{n+1}-ru_{i-1}^{n+1}=ru_{i+1}^{n}+(1-2r)u_{i}^{n}+ru_{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 Verfahrens 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.

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: der bekannten Diffusionsgleichung (gewählt als Konstante $D_{x}$), einer advektiven Komponente (das System entwickelt sich infolge eines Vektorfeldes im Raum), welche als Konstante $U_{x}$ gewählt wird, und einer lateralen Ü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 bzw. 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}-2C_{i}^{j+1}+C_{i-1}^{j+1})+(C_{i+1}^{j}-2C_{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})$

Die folgenden Konstanten werden 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}}\quad {\text{und}}\quad \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}}\quad {\text{und}}\quad \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 4×4-Nullmatrix entsprechen. Es ist zu beachten, dass AA und BB jeweils die Größe 12×12 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}}\quad {\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 12×1-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{aligned}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{aligned}}$

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 $\mu$ 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{aligned}&(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{aligned}}$

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. 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.

Weblinks

• Modul für parabolische partielle DGL (engl.) (Memento vom 27. September 2013 im Internet Archive)