Crank–Nicolson method

From Wikipedia, the free encyclopedia

In the mathematical subfield numerical analysis, the Crank–Nicolson method is a finite difference method used for numerically solving the heat equation and similar partial differential equations.[1] It is a second-order method in time, implicit in time, and is numerically stable. The method was developed by John Crank and Phyllis Nicolson in the mid 20th century.

For diffusion equations (and many other equations), it can be shown the Crank–Nicolson method is unconditionally stable. However, the approximate solutions can still contain (decaying) spurious oscillations if the ratio of time step to the square of space step is large (typically larger than 1 / 2). For this reason, whenever large time steps or high spatial resolution is necessary, the less accurate backward Euler method is often used, which is both stable and immune to oscillations.

Contents

[edit] The method

The Crank–Nicolson stencil on a 1D problem.
The Crank–Nicolson stencil on a 1D problem.

The Crank–Nicolson method is based on central difference in space, and the trapezoidal rule in time, giving second-order convergence in time. Equivalently, it is the average of forward Euler and backward Euler in time. For example, in one dimension, if the partial differential equation is

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

then, letting u(i \Delta x, n \Delta t) = u_{i}^{n}\,, the Crank–Nicolson method is the average of the forward Euler method at n and the backward Euler method at 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{(forward Euler)}
\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{(backward Euler)}
\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)}

The function F must be discretized spatially with a central difference.

Note that this is an implicit method: to get the "next" value of u in time, a system of algebraic equations must be solved. If the partial differential equation is nonlinear, the discretization will also be nonlinear so that advancing in time will involve the solution of a system of nonlinear algebraic equations, though linearizations are possible. In many problems, especially linear diffusion, the algebraic problem is tridiagonal and may be efficiently solved with the tridiagonal matrix algorithm, avoiding a costly full matrix inversion.

[edit] Example: 1D diffusion

The Crank–Nicolson method is often applied to diffusion problems. As an example, for linear diffusion,

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

whose Crank–Nicolson discretization is then:

\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)

or, letting 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}\,

which is a tridiagonal problem, so that u_{i}^{n + 1}\, may be efficiently solved for using the tridiagonal matrix algorithm in favor of a much more costly matrix inversion.

A quasilinear equation, such as (this is a minimalistic example and not general)

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

would lead to a nonlinear system of algebraic equations which could not be easily solved as above; however, it is possible in some cases to linearize the problem by using the old value for a, that is a_{i}^{n}(u)\, instead of a_{i}^{n + 1}(u)\,. Other times, it may be possible to estimate a_{i}^{n + 1}(u)\, using an explicit method and maintain stability.

[edit] Example: 1D diffusion with advection for steady flow, with multiple channel connections

This is a solution usually employed for many purposes when there's a contamination problem in streams or rivers under steady flow conditions but the only information is given in one dimension. Sometimes is needed some kind of info or data about the behavior in the cross section so this is one of many ways to get some extra data without getting into a model of two or three dimensions.

What is modeled here is the concentration of a solute contaminant in water for example. This problem is composed by three parts the known diffusion equation (Dx chosen as constant), plus an advective component which means the system is evolving in space due to a velocity field (Ux which in this case we choose as constant). The third part of the equation is due to a lateral interaction between longitudinal channels (k).

<0>\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)

where C is the concentration of the contaminant subscripts N & M correspond to previous and next channel.

The Crank–Nicolson scheme of solution (knowing i as position in space & j as time interval) comes like this:

<1> \frac{\partial C}{\partial t} = \frac{C_{i}^{j + 1} - C_{i}^{j}}{\Delta t}
<2>\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>\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> C= \frac{1}{2} (C_{i}^{j+1} + C_{i}^{j})
<5> C_N= \frac{1}{2} (C_{Ni}^{j+1} + C_{Ni}^{j})
<6> C_M= \frac{1}{2} (C_{Mi}^{j+1} + C_{Mi}^{j})


Now we create the following constants that will be helpful doing the algebra:

\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}


replacing <1>,<2>,<3>,<4>,<5>,<6>, α, β & λ into <0> and putting the terms for the new time in the left (j+1) and the term for the present time in the right (j) we get this expression:

-\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 the case we are trying to model the first channel, it can only be in contact with the following channel (M), so the expression is simplified by something like this:


-(\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}


in the same way if the model is for the last channel, it can only be in contact with the previous channel (N), so the expression is simplified by something like this:


-\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}


in order to solve this kind linear of system of equations we must now see that boundary conditions must be given first to the beginning of the channels:


C_0^{j}: initial condition for the channel at present time step
C_{0}^{j+1}: initial condition for the channel at next time step
C_{N0}^{j}: initial condition for the previous channel to the one analized at present time step
C_{M0}^{j}: initial condition for the next channel to the one analized at present time step


and for the last cell of the channels (z) we must see that the most convenient condition becomes an adiabatic one which means


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


condition satisfied if and only if (regardless of a null value)

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


Now let's see what happens if the problem is solved (in a matrix form) using 3 channels and 5 nodes (including the initial boundary condition) we can express this as a linear system problem like this one:


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

where

\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}   and   \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}


now we must realize that AA and BB should be arrays made of 4 different subarrays (remember that only three channels are considered for this example but it covers the main part discussed above).

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


where the elements mentioned above correspond to the next arrays and an additional 4x4 full of zeros. Must be noted that the AA and BB arrays had size 12x12:


\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}   &  
\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 d vector is just for the assignment of the boundary conditions. In this example is a 12x1 vector:

\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}

One must the iterate the following equation to solve the system for any time

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

[edit] Example: 2D diffusion

When extending into two dimensions on a uniform Cartesian grid, the derivation is similar and the results may lead to a system of band-diagonal equations rather than tridiagonal ones. The two-dimensional heat equation

\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 of

\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}

assuming that a square grid is used so that Δx = Δy. This equation can be simplified somewhat by rearranging terms and using the CFL number

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

For the Crank–Nicolson numerical scheme, a low CFL number is not required for stability, however it is required for numerical accuracy. We can now write the scheme as:

\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}

[edit] Application in financial mathematics

Because a number of other phenomena can be modeled with the heat equation (often called the diffusion equation in financial mathematics), the Crank–Nicolson method has been applied to those areas as well. Particularly, the Black-Scholes option pricing model's differential equation can be transformed into the heat equation, and thus option pricing numerical solutions can be obtained with the Crank–Nicolson method. The importance of that comes from the extensions of the option pricing model that are not able to be represented with a closed form analytic solution; they can still offer numerical solutions. However, for non-smooth final conditions (which happen for most financial instruments), the Crank–Nicolson method is not satisfactory as numerical oscillations are not damped. For vanilla options, this results in oscillation in the gamma value around the strike price. Therefore, special damping initialization steps are necessary (e.g., fully implicit finite difference method).

[edit] See also

[edit] References

  1. ^ Tuncer Cebeci (2002). Convective Heat Transfer. Springer. ISBN 0966846141. 
  • Crank J. and Nicolson P. (1947) "A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type". Proceedings of the Cambridge Philosophical Society 43, 50–64.
  • Wilmott P., Howison S., Dewynne J. (1995) The Mathematics of Financial Derivatives:A Student Introduction. Cambridge University Press.

[edit] External links

Languages