Crank-Nicolson method

In numerical mathematics, the Crank-Nicolson method is a finite difference method for solving the heat conduction equation and similar partial differential equations . It is an implicit procedure of the 2nd order and numerically stable . The process was developed by John Crank and Phyllis Nicolson in the middle of the 20th century .

For the heat conduction equation and many other equations it can be shown that the Crank-Nicolson method is numerically stable without conditions. Nevertheless, the approximated solutions can contain disruptive oscillations if the quotient of the time difference and the distance square is large (typically greater than ). In this case, the less precise Euler backward method is often used, which is numerically stable and insensitive to interference. ${\displaystyle 1/2}$

The procedure

Difference star ( stencil ) for the Crank – Nicolson method in the 1D problem. j is the place and n is the time.

The Crank-Nicolson method is based on the trapezoid rule in time. This is equivalent to the mean value of the Euler forward method and the Euler backward method .

Is the partial differential equation in the one-dimensional case

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

given and denote , then the Crank-Nicolson method is the mean value of the Euler forward method for value and the Euler backward method for : ${\displaystyle u(i\Delta x,n\Delta t)=u_{i}^{n}\,}$ ${\displaystyle n}$ ${\displaystyle n+1}$

${\displaystyle {\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)}}}$

${\displaystyle {\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)}}}$

${\displaystyle {\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)}}}$

The indices and stand for the place or time. The function must be spatially discretized, e.g. B. with the finite difference method , the finite volume method or the finite element method . ${\displaystyle i}$ ${\displaystyle n}$ ${\displaystyle F}$

It should be noted that this is an implicit method. In order to get the temporally "closest" value of , an algebraic system of equations has to be solved. If the partial differential equation is not linear, the discretization will also not be linear, so that a non-linear algebraic system of equations must be solved, although linearizations are possible. With many problems, especially with linear diffusion, the algebraic problem is tridiagonal and is best solved with the fast Thomas algorithm in order to avoid a complex inversion of the matrix. ${\displaystyle u}$

Example: One-dimensional diffusion

The Crank-Nicolson method is often used to solve diffusion problems. An example of linear diffusion is

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

Using the finite difference method for the spatial discretization of the right side, the Crank – Nicolson discretization is then:

${\displaystyle {\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)}$

or, be : ${\displaystyle r={\frac {a\Delta t}{2(\Delta x)^{2}}}}$

${\displaystyle -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}\,}$

this is a tridiagonal problem, so that it should be solved with the Thomas algorithm in order to avoid a complex inversion of the resulting tridiagonal Toeplitz matrix . ${\displaystyle u_{i}^{n+1}\,}$

A quasi-linear equation like (this is a simplified example and not generally applicable)

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

would lead to a non-linear system of algebraic equations that are not as easy to solve as above. However, in some cases it is possible to linearize the problem using the old value of ; this procedure is also referred to as "lagging the coefficients" . This is taking place . In other cases it may be possible to estimate using an explicit method and maintaining stability. ${\displaystyle a}$ ${\displaystyle a_{i}^{n}(u)\,}$ ${\displaystyle a_{i}^{n+1}(u)\,}$ ${\displaystyle a_{i}^{n+1}(u)\,}$

Alternatively, iterative updating of the coefficients enables the order of the method to be maintained. Here, "Lagging the Coefficients" applied in a first step to solve the system of equations for the next time step. With the new temperatures now known, the coefficients are evaluated and the solution is calculated again. This procedure is continued until the difference between the solutions of two successive iterations is small.

Example: One-dimensional diffusion with advection for stationary flows

This solution is usually used when there is contamination in rivers with steady currents, but information is only given in one dimension. Often times the problem can be simplified into a one-dimensional problem and still yield useful information.

In the following, a dissolved foreign body is modeled in water. This problem is composed of three parts: the well-known diffusion equation (chosen as a constant ), an advective component (the system develops as a result of a vector field in space), which is chosen as a constant , and a lateral superposition between longitudinal channels ( ) ${\displaystyle D_{x}}$ ${\displaystyle U_{x}}$ ${\displaystyle k}$

${\displaystyle {\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 (*)}$

where the concentration of the foreign body is in the water and the indices and correspond to the previous and the next channel. ${\displaystyle C}$ ${\displaystyle N}$ ${\displaystyle M}$

The Crank-Nicolson method (where stands for place and time) transforms each component of the partial differential equation as follows: ${\displaystyle i}$ ${\displaystyle j}$

${\displaystyle 1.\quad {\frac {\partial C}{\partial t}}={\frac {C_{i}^{j+1}-C_{i}^{j}}{\Delta t}}}$

${\displaystyle 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)}$

${\displaystyle 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)}$

${\displaystyle 4.\quad C={\frac {1}{2}}(C_{i}^{j+1}+C_{i}^{j})}$

${\displaystyle 5.\quad C_{N}={\frac {1}{2}}(C_{Ni}^{j+1}+C_{Ni}^{j})}$

${\displaystyle 6.\quad C_{M}={\frac {1}{2}}(C_{Mi}^{j+1}+C_{Mi}^{j})}$

The following constants are defined to simplify the calculation:

${\displaystyle \lambda ={\frac {D_{x}\Delta t}{2\Delta x^{2}}}}$

${\displaystyle \alpha ={\frac {U_{x}\Delta t}{4\Delta x}}}$

${\displaystyle \beta ={\frac {k\Delta t}{2}}}$

Now 1. to 6. ,, and are inserted in (*). Then the “future” terms ( ) are moved to the left and the “current” ( ) to the right. It turns out ${\displaystyle \alpha }$ ${\displaystyle \beta }$ ${\displaystyle \lambda }$ ${\displaystyle j+1}$ ${\displaystyle j}$

${\displaystyle -\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}.}$

In order to model the first channel, it is stated that it can only be in contact with the following channel ( M ), so that the expression is simplified to

${\displaystyle -(\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}.}$

The last channel is modeled in the same way. It is stated that it can only be in contact with the previous channel ( N ), so the expression becomes too simplified

${\displaystyle -\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}.}$

To solve this system of linear equations, boundary conditions must be given at the beginning of the channels:

${\displaystyle C_{0}^{j}}$: Initial value for the current channel

${\displaystyle C_{0}^{j+1}}$: Initial value for the channel of the next time interval

${\displaystyle C_{N0}^{j}}$: Initial value for the previous channel of the current time interval

${\displaystyle C_{M0}^{j}}$: Start value for the following channel of the current time interval

For the last cell of the channels ( ) the most favorable condition becomes adiabatic, that is ${\displaystyle z}$

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

This condition is fulfilled if (and only then (apart from the value 0))

${\displaystyle C_{i+1}^{j+1}=C_{i-1}^{j+1}.}$

We now want to solve the problem (in matrix form) for the case with three channels and five nodes (including the initial value conditions). Expressed as a linear problem:

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

in which

${\displaystyle \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}}.}$

It is stated that AA and BB should be fields with four different field components (as a reminder, only three channels have been considered for this example, but it still covers most of the above):

${\displaystyle \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}},}$  

where the elements mentioned above correspond to the next fields and an additional 4 × 4 zero matrix . It should be noted that AA and BB are each 12 × 12 in size:

${\displaystyle \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}},}$

${\displaystyle \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}},}$

${\displaystyle \mathbf {AA3} ={\begin{bmatrix}-\beta &0&0&0\\0&-\beta &0&0\\0&0&-\beta &0\\0&0&0&-\beta \end{bmatrix}},}$

${\displaystyle \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}}}$

${\displaystyle \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}}.}$

The vector is used to keep the boundary conditions. In this example it is a 12 × 1 vector: ${\displaystyle d}$

${\displaystyle \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}}.}$

To find the concentration at any point in time, the following equation must be solved iteratively:

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

Example: two-dimensional diffusion

If two dimensions are considered, the derivation is analogous and the results provide a system of band-diagonal rather than tridiagonal equations. The two-dimensional heat conduction equation

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

can be solved with the Crank-Nicolson discretization

${\displaystyle {\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}}}$

assumed that a rectangular coordinate network is used, so that . This equation can be simplified a little by transforming the term and using the CFL number , ${\displaystyle \Delta x=\Delta y}$ ${\displaystyle \mu }$

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

A low CFL number is not necessary for the stability of the Crank-Nicolson numerical scheme, but it is necessary for numerical accuracy. The scheme can now be represented as

${\displaystyle {\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}}}$

Application in financial mathematics

Because a number of other phenomena can also be modeled with the help of the heat conduction equation ( often called diffusion equation in financial mathematics ), the Crank-Nicolson method is also used in these areas. In particular, the differential equations of the Black-Scholes model can be converted into the diffusion equation and solved using the Crank-Nicolson method. The importance of this comes from the spread of the option price model, which cannot be represented with a self-contained analytical solution; numerical solutions can still offer themselves. However, the Crank-Nicolson method is unsatisfactory in uneven financial conditions (which is common to most financial instruments) because numerical oscillations are not dampened. Therefore special damping initializations are necessary.

Individual evidence

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 . In: Springer (Ed.): Advances in Computational Mathematics . 6th vol., No. 1, December 1996, ISSN  1019-7168 , pp. 207-226. doi : 10.1007 / BF02127704 .
3. ↑ Pletcher, Richard H., John C. Tannehill, Dale Anderson: Computational fluid mechanics and heat transfer. CRC Press, 2012.
4. ↑ Wilmott, P .; Howison, S .; Dewynne, J .: The Mathematics of Financial Derivatives: A Student Introduction . Cambridge Univ. Press, 1995, ISBN 0521497892 .

Web links

• Module (Engl.) For parabolic partial differential equation ( Memento of 27 September 2013 Internet Archive )