Menu
Home
People
Places
Arts
History
Plants & Animals
Science
Life & Culture
Technology
Reference.org
Crank–Nicolson method
open-in-new
<h2 id="principle">Principle</h2> <p>The Crank–Nicolson method is based on the <a href="/facts/Trapezoidal_rule_(differential_equations)/vlUlTgYv">trapezoidal rule</a>, giving second-order convergence in time. For linear equations, the trapezoidal rule is equivalent to the <a href="/facts/Midpoint_method/7So6u6Yo">implicit midpoint method</a>—the simplest example of a <a href="/facts/Gauss%25E2%2580%2593Legendre_method/o6gsDLuT">Gauss–Legendre</a> implicit Runge–Kutta method—which also has the property of being a <a href="/facts/Geometric_integrator/Yz76jzPK">geometric integrator</a>. For example, in one dimension, suppose the <a href="/facts/Partial_differential_equation/DyPsBlv7">partial differential equation</a> is ∂ u ∂ t = F ( u , x , t , ∂ u ∂ x , ∂ 2 u ∂ x 2 ) . {\displaystyle {\frac {\partial u}{\partial t}}=F\left(u,x,t,{\frac {\partial u}{\partial x}},{\frac {\partial ^{2}u}{\partial x^{2}}}\right).} </p><p>Letting u ( i Δ x , n Δ t ) = u i n {\displaystyle u(i\Delta x,n\Delta t)=u_{i}^{n}} and F i n = F {\displaystyle F_{i}^{n}=F} evaluated for i , n {\displaystyle i,n} and u i n {\displaystyle u_{i}^{n}} , the equation for Crank–Nicolson method is a combination of the <a href="/facts/Forward_Euler_method/7yF5B4z9">forward Euler method</a> at n {\displaystyle n} and the <a href="/facts/Backward_Euler_method/lccBRRYY">backward Euler method</a> at n + 1 {\displaystyle n+1} (note, however, that the method itself is <i>not</i> simply the average of those two methods, as the backward Euler equation has an implicit dependence on the solution): </p> <table><tbody><tr><td> u i n + 1 − u i n Δ t = F i n ( u , x , t , ∂ u ∂ x , ∂ 2 u ∂ x 2 ) {\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)} </td><td>forward Euler</td></tr><tr><td> u i n + 1 − u i n Δ t = F i n + 1 ( u , x , t , ∂ u ∂ x , ∂ 2 u ∂ x 2 ) {\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)} </td><td>backward Euler</td></tr><tr><td> u i n + 1 − u i n Δ t = 1 2 [ F i n + 1 ( u , x , t , ∂ u ∂ x , ∂ 2 u ∂ x 2 ) + F i n ( u , x , t , ∂ u ∂ x , ∂ 2 u ∂ x 2 ) ] {\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]} </td><td>Crank–Nicolson</td></tr></tbody></table> <p>Note that this is an <i>implicit method</i>: to get the "next" value of u {\displaystyle u} in time, a system of algebraic equations must be solved. If the partial differential equation is nonlinear, the <a href="/facts/Temporal_discretization/nXPjUf10">discretization</a> 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 <a href="/facts/Tridiagonal/1SK5A0rP">tridiagonal</a> and may be efficiently solved with the <a href="/facts/Tridiagonal_matrix_algorithm/NREF5x1c">tridiagonal matrix algorithm</a>, which gives a fast O ( N ) {\displaystyle {\mathcal {O}}(N)} direct solution, as opposed to the usual O ( N 3 ) {\displaystyle {\mathcal {O}}(N^{3})} for a full matrix, in which N {\displaystyle N} indicates the matrix size. </p> <h2 id="example-1d-diffusion">Example: 1D diffusion</h2> <p>The Crank–Nicolson method is often applied to <a href="/facts/Diffusion_equation/l3CUuLqF">diffusion problems</a>. As an example, for linear diffusion, </p> ∂ u ∂ t = a ∂ 2 u ∂ x 2 , {\displaystyle {\frac {\partial u}{\partial t}}=a{\frac {\partial ^{2}u}{\partial x^{2}}},} <p>applying a <a href="/facts/Finite_difference/aY1txFhj">finite difference</a> spatial discretization for the right-hand side, the Crank–Nicolson discretization is then </p> u i n + 1 − u i n Δ t = a 2 ( Δ x ) 2 ( ( 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 ) ) {\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)} <p>or, letting r = a Δ t 2 ( Δ x ) 2 {\displaystyle r={\frac {a\Delta t}{2(\Delta x)^{2}}}} , </p> − 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 . {\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}.} <p>Given that the terms on the right-hand side of the equation are known, this is a <a href="/facts/Tridiagonal/1SK5A0rP">tridiagonal</a> problem, so that u i n + 1 {\displaystyle u_{i}^{n+1}} may be efficiently solved by using the <a href="/facts/Tridiagonal_matrix_algorithm/NREF5x1c">tridiagonal matrix algorithm</a> in favor over the much more costly <a href="/facts/Matrix_inversion/fPqXk3V8">matrix inversion</a>. </p><p>A quasilinear equation, such as (this is a minimalistic example and not general) </p> ∂ u ∂ t = a ( u ) ∂ 2 u ∂ x 2 , {\displaystyle {\frac {\partial u}{\partial t}}=a(u){\frac {\partial ^{2}u}{\partial x^{2}}},} <p>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 {\displaystyle a} , that is, a i n ( u ) {\displaystyle a_{i}^{n}(u)} instead of a i n + 1 ( u ) {\displaystyle a_{i}^{n+1}(u)} . Other times, it may be possible to estimate a i n + 1 ( u ) {\displaystyle a_{i}^{n+1}(u)} using an explicit method and maintain stability. </p> <h2 id="example-1d-diffusion-with-advection-for-steady-flow-with-multiple-channel-connections">Example: 1D diffusion with advection for steady flow, with multiple channel connections</h2> <p>This is a solution usually employed for many purposes when there is a contamination problem in streams or rivers under steady flow conditions, but information is given in one dimension only. Often the problem can be simplified into a 1-dimensional problem and still yield useful information. </p><p>Here we model the concentration of a solute contaminant in water. This problem is composed of three parts: the known diffusion equation ( D x {\displaystyle D_{x}} chosen as constant), an <a href="/facts/Advection/3rn50V2r">advective</a> component (which means that the system is evolving in space due to a velocity field), which we choose to be a constant U x {\displaystyle U_{x}} , and a lateral interaction between longitudinal channels ( k {\displaystyle k} ): </p> <table><tbody><tr><td> ∂ C ∂ t = D x ∂ 2 C ∂ x 2 − U x ∂ C ∂ x − k ( C − C N ) − k ( C − C M ) , {\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}),} </td> <td></td> <td>1</td></tr></tbody></table> <p>where C {\displaystyle C} is the concentration of the contaminant, and subscripts N {\displaystyle N} and M {\displaystyle M} correspond to <i>previous</i> and <i>next</i> channel. </p><p>The Crank–Nicolson method (where i {\displaystyle i} represents position, and j {\displaystyle j} time) transforms each component of the PDE into the following: </p> <table><tbody><tr><td> ∂ C ∂ t ⇒ C i j + 1 − C i j Δ t , {\displaystyle {\frac {\partial C}{\partial t}}\Rightarrow {\frac {C_{i}^{j+1}-C_{i}^{j}}{\Delta t}},} </td> <td></td> <td>2</td></tr></tbody></table> <table><tbody><tr><td> ∂ 2 C ∂ x 2 ⇒ 1 2 ( Δ x ) 2 ( ( 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 ) ) , {\displaystyle {\frac {\partial ^{2}C}{\partial x^{2}}}\Rightarrow {\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),} </td> <td></td> <td>3</td></tr></tbody></table> <table><tbody><tr><td> ∂ C ∂ x ⇒ 1 2 ( ( C i + 1 j + 1 − C i − 1 j + 1 ) 2 ( Δ x ) + ( C i + 1 j − C i − 1 j ) 2 ( Δ x ) ) , {\displaystyle {\frac {\partial C}{\partial x}}\Rightarrow {\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),} </td> <td></td> <td>4</td></tr></tbody></table> <table><tbody><tr><td> C ⇒ 1 2 ( C i j + 1 + C i j ) , {\displaystyle C\Rightarrow {\frac {1}{2}}(C_{i}^{j+1}+C_{i}^{j}),} </td> <td></td> <td>5</td></tr></tbody></table> <table><tbody><tr><td> C N ⇒ 1 2 ( C N i j + 1 + C N i j ) , {\displaystyle C_{N}\Rightarrow {\frac {1}{2}}(C_{Ni}^{j+1}+C_{Ni}^{j}),} </td> <td></td> <td>6</td></tr></tbody></table> <table><tbody><tr><td> C M ⇒ 1 2 ( C M i j + 1 + C M i j ) . {\displaystyle C_{M}\Rightarrow {\frac {1}{2}}(C_{Mi}^{j+1}+C_{Mi}^{j}).} </td> <td></td> <td>7</td></tr></tbody></table> <p>Now we create the following constants to simplify the algebra: </p> λ = D x Δ t 2 Δ x 2 , {\displaystyle \lambda ={\frac {D_{x}\,\Delta t}{2\,\Delta x^{2}}},} α = U x Δ t 4 Δ x , {\displaystyle \alpha ={\frac {U_{x}\,\Delta t}{4\,\Delta x}},} β = k Δ t 2 , {\displaystyle \beta ={\frac {k\,\Delta t}{2}},} <p>and substitute (2), (3), (4), (5), (6), (7), α {\displaystyle \alpha } , β {\displaystyle \beta } and λ {\displaystyle \lambda } into (1). We then put the <i>new time</i> terms on the left ( j + 1 {\displaystyle j+1} ) and the <i>present time</i> terms on the right ( j {\displaystyle j} ) to get </p> − β C N i j + 1 − ( λ + α ) C i − 1 j + 1 + ( 1 + 2 λ + 2 β ) C i j + 1 − ( λ − α ) C i + 1 j + 1 − β C M i j + 1 = {\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}={}} β C N i j + ( λ + α ) C i − 1 j + ( 1 − 2 λ − 2 β ) C i j + ( λ − α ) C i + 1 j + β C M i j . {\displaystyle \qquad \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}.} <p>To model the <i>first</i> channel, we realize that it can only be in contact with the following channel ( M {\displaystyle M} ), so the expression is simplified to </p> − ( λ + α ) C i − 1 j + 1 + ( 1 + 2 λ + β ) C i j + 1 − ( λ − α ) C i + 1 j + 1 − β C M i j + 1 = {\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}={}} + ( λ + α ) C i − 1 j + ( 1 − 2 λ − β ) C i j + ( λ − α ) C i + 1 j + β C M i j . {\displaystyle \qquad {}+(\lambda +\alpha )C_{i-1}^{j}+(1-2\lambda -\beta )C_{i}^{j}+(\lambda -\alpha )C_{i+1}^{j}+\beta C_{Mi}^{j}.} <p>In the same way, to model the <i>last</i> channel, we realize that it can only be in contact with the previous channel ( N {\displaystyle N} ), so the expression is simplified to </p> − β C N i j + 1 − ( λ + α ) C i − 1 j + 1 + ( 1 + 2 λ + β ) C i j + 1 − ( λ − α ) C i + 1 j + 1 = {\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}={}} β C N i j + ( λ + α ) C i − 1 j + ( 1 − 2 λ − β ) C i j + ( λ − α ) C i + 1 j . {\displaystyle \qquad \beta C_{Ni}^{j}+(\lambda +\alpha )C_{i-1}^{j}+(1-2\lambda -\beta )C_{i}^{j}+(\lambda -\alpha )C_{i+1}^{j}.} <p>To solve this linear system of equations, we must now see that boundary conditions must be given first to the beginning of the channels: </p> C 0 j {\displaystyle C_{0}^{j}} : boundary condition for the channel at present time step, C 0 j + 1 {\displaystyle C_{0}^{j+1}} : boundary condition for the channel at next time step, C N 0 j {\displaystyle C_{N0}^{j}} : boundary condition for the previous channel to the one analyzed at present time step, C M 0 j {\displaystyle C_{M0}^{j}} : boundary condition for the next channel to the one analyzed at present time step. <p>For the last cell of the channels ( z {\displaystyle z} ), the most convenient condition becomes an adiabatic one, so </p> ∂ C ∂ x | x = z = C i + 1 − C i − 1 2 Δ x = 0. {\displaystyle \left.{\frac {\partial C}{\partial x}}\right|_{x=z}={\frac {C_{i+1}-C_{i-1}}{2\,\Delta x}}=0.} <p>This condition is satisfied if and only if (regardless of a null value) </p> C i + 1 j + 1 = C i − 1 j + 1 . {\displaystyle C_{i+1}^{j+1}=C_{i-1}^{j+1}.} <p>Let us solve this problem (in a matrix form) for the case of 3 channels and 5 nodes (including the initial boundary condition). We express this as a linear system problem: </p> A A C j + 1 = B B C j + d , {\displaystyle \mathbf {AA} \,\mathbf {C^{j+1}} =\mathbf {BB} \,\mathbf {C^{j}} +\mathbf {d} ,} <p>where </p> C j + 1 = [ 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 ] , C j = [ 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 ] . {\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 \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}}.} <p>Now we must realize that AA and BB should be arrays made of four different subarrays (remember that only three channels are considered for this example, but it covers the main part discussed above): </p> A A = [ A A 1 A A 3 0 A A 3 A A 2 A A 3 0 A A 3 A A 1 ] , B B = [ B B 1 − A A 3 0 − A A 3 B B 2 − A A 3 0 − A A 3 B B 1 ] , {\displaystyle \mathbf {AA} ={\begin{bmatrix}AA1&AA3&0\\AA3&AA2&AA3\\0&AA3&AA1\end{bmatrix}},\quad \mathbf {BB} ={\begin{bmatrix}BB1&-AA3&0\\-AA3&BB2&-AA3\\0&-AA3&BB1\end{bmatrix}},} <p>where the elements mentioned above correspond to the next arrays, and an additional 4×4 full of zeros. Please note that the sizes of AA and BB are 12×12: </p> A A 1 = [ ( 1 + 2 λ + β ) − ( λ − α ) 0 0 − ( λ + α ) ( 1 + 2 λ + β ) − ( λ − α ) 0 0 − ( λ + α ) ( 1 + 2 λ + β ) − ( λ − α ) 0 0 − 2 λ ( 1 + 2 λ + β ) ] , {\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}},} A A 2 = [ ( 1 + 2 λ + 2 β ) − ( λ − α ) 0 0 − ( λ + α ) ( 1 + 2 λ + 2 β ) − ( λ − α ) 0 0 − ( λ + α ) ( 1 + 2 λ + 2 β ) − ( λ − α ) 0 0 − 2 λ ( 1 + 2 λ + 2 β ) ] , {\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}},} A A 3 = [ − β 0 0 0 0 − β 0 0 0 0 − β 0 0 0 0 − β ] , {\displaystyle \mathbf {AA3} ={\begin{bmatrix}-\beta &0&0&0\\0&-\beta &0&0\\0&0&-\beta &0\\0&0&0&-\beta \end{bmatrix}},} B B 1 = [ ( 1 − 2 λ − β ) ( λ − α ) 0 0 ( λ + α ) ( 1 − 2 λ − β ) ( λ − α ) 0 0 ( λ + α ) ( 1 − 2 λ − β ) ( λ − α ) 0 0 2 λ ( 1 − 2 λ − β ) ] , {\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}},} B B 2 = [ ( 1 − 2 λ − 2 β ) ( λ − α ) 0 0 ( λ + α ) ( 1 − 2 λ − 2 β ) ( λ − α ) 0 0 ( λ + α ) ( 1 − 2 λ − 2 β ) ( λ − α ) 0 0 2 λ ( 1 − 2 λ − 2 β ) ] . {\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}}.} <p>The d vector here is used to hold the boundary conditions. In this example it is a 12×1 vector: </p> d = [ ( λ + α ) ( C 10 j + 1 + C 10 j ) 0 0 0 ( λ + α ) ( C 20 j + 1 + C 20 j ) 0 0 0 ( λ + α ) ( C 30 j + 1 + C 30 j ) 0 0 0 ] . {\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}}.} <p>To find the concentration at any time, one must iterate the following equation: </p> C j + 1 = A A − 1 ( B B C j + d ) . {\displaystyle \mathbf {C^{j+1}} =\mathbf {AA} ^{-1}(\mathbf {BB} \,\mathbf {C^{j}} +\mathbf {d} ).} <h2 id="example-2d-diffusion">Example: 2D diffusion</h2> <p>When extending into two dimensions on a uniform <a href="/facts/Cartesian_grid/8vDdL4UB">Cartesian grid</a>, the derivation is similar and the results may lead to a system of <a href="/facts/Banded_matrix/Iz803N35">band-diagonal</a> equations rather than <a href="/facts/Tridiagonal_matrix/1SK5A0rP">tridiagonal</a> ones. The two-dimensional heat equation </p> ∂ u ∂ t = a ∇ 2 u , {\displaystyle {\frac {\partial u}{\partial t}}=a\,\nabla ^{2}u,} ∂ u ∂ t = a ( ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 ) {\displaystyle {\frac {\partial u}{\partial t}}=a\left({\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}\right)} <p>can be solved with the Crank–Nicolson discretization of </p> u i , j n + 1 = u i , j n + 1 2 a Δ t ( Δ x ) 2 [ ( u i + 1 , j n + 1 + u i − 1 , j n + 1 + u i , j + 1 n + 1 + u i , j − 1 n + 1 − 4 u i , j n + 1 ) + ( u i + 1 , j n + u i − 1 , j n + u i , j + 1 n + u i , j − 1 n − 4 u i , j n ) ] , {\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})\\&+(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}}} <p>assuming that a square grid is used, so that Δ x = Δ y {\displaystyle \Delta x=\Delta y} . This equation can be simplified somewhat by rearranging terms and using the <a href="/facts/Courant_number/f4r00Wca">CFL number</a> </p> μ = a Δ t ( Δ x ) 2 . {\displaystyle \mu ={\frac {a\,\Delta t}{(\Delta x)^{2}}}.} <p>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 </p> ( 1 + 2 μ ) u i , j n + 1 − μ 2 ( u i + 1 , j n + 1 + u i − 1 , j n + 1 + u i , j + 1 n + 1 + u i , j − 1 n + 1 ) {\displaystyle (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)} = ( 1 − 2 μ ) u i , j n + μ 2 ( u i + 1 , j n + u i − 1 , j n + u i , j + 1 n + u i , j − 1 n ) . {\displaystyle \qquad =(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).} <p>Solving such a linear system is costly. Hence an <a href="/facts/Alternating-direction_implicit_method/uV9e0ccz">alternating-direction implicit method</a> can be implemented to solve the numerical PDE, whereby one dimension is treated implicitly, and other dimension explicitly for half of the assigned time step and conversely for the remainder half of the time step. The benefit of this strategy is that the implicit solver only requires a <a href="/facts/Tridiagonal_matrix_algorithm/NREF5x1c">tridiagonal matrix algorithm</a> to be solved. The difference between the true Crank–Nicolson solution and ADI approximated solution has an order of accuracy of O ( Δ t 4 ) {\displaystyle O(\Delta t^{4})} and hence can be ignored with a sufficiently small time step.<a class="footnote-ref" id="fnref:4" href="#fn:4"><sup>4</sup></a> </p> <h2 id="cranknicolson-for-nonlinear-problems">Crank–Nicolson for nonlinear problems</h2> <p>Because the Crank–Nicolson method is <a href="/facts/Implicit_function/UJmQ2Cyk">implicit</a>, it is generally impossible to solve exactly. Instead, an iterative technique should be used to converge to the solution. One option is to use <a href="/facts/Newton%2527s_method/VlRhI2FL">Newton's method</a> to converge on the prediction, but this requires the computation of the <a href="/facts/Jacobian_matrix_and_determinant/0tNv6o8j">Jacobian</a>. For a high-dimensional system like those in <a href="/facts/Computational_fluid_dynamics/zuWdtNJp">computational fluid dynamics</a> or <a href="/facts/Numerical_relativity/R7D2PEqT">numerical relativity</a>, it may be infeasible to compute this Jacobian. </p><p>A Jacobian-free alternative is <a href="/facts/Fixed-point_iteration/05dQGazd">fixed-point iteration</a>. If f {\displaystyle f} is the velocity of the system, then the Crank–Nicolson prediction will be a fixed point of the map Φ ( x ) = x 0 + h 2 [ f ( x 0 ) + f ( x ) ] . {\displaystyle \Phi (x)=x_{0}+{\frac {h}{2}}\left[f(x_{0})+f(x)\right].} If the map iteration x ( i + 1 ) = Φ ( x ( i ) ) {\displaystyle x^{(i+1)}=\Phi (x^{(i)})} does not converge, the parameterized map Θ ( x , α ) = α x + ( 1 − α ) Φ ( x ) {\displaystyle \Theta (x,\alpha )=\alpha x+(1-\alpha )\Phi (x)} , with α ∈ ( 0 , 1 ) {\displaystyle \alpha \in (0,1)} , may be better behaved. In expanded form, the update formula is </p> x i + 1 = α x i + ( 1 − α ) [ x 0 + h 2 ( f ( x 0 ) + f ( x i ) ) ] , {\displaystyle x^{i+1}=\alpha x^{i}+(1-\alpha )\left[x_{0}+{\frac {h}{2}}\left(f(x_{0})+f(x^{i})\right)\right],} <p>where x i {\displaystyle x^{i}} is the current guess and x i − 1 {\displaystyle x_{i-1}} is the previous time-step. </p><p>Even for high-dimensional systems, iteration of this map can converge surprisingly quickly. </p> <h2 id="application-in-financial-mathematics">Application in financial mathematics</h2> <p class="note">Further information: <a href="/facts/Finite_difference_methods_for_option_pricing/kLGugeug">Finite difference methods for option pricing</a></p> <p>Because a number of other phenomena can be <a href="/facts/Mathematical_model/CGhyuxgz">modeled</a> with the <a href="/facts/Heat_equation/2JPPr5GG">heat equation</a> (often called the diffusion equation in <a href="/facts/Financial_mathematics/P9otEltQ">financial mathematics</a>), the Crank–Nicolson method has been applied to those areas as well.<a class="footnote-ref" id="fnref:5" href="#fn:5"><sup>5</sup></a> Particularly, the <a href="/facts/Black%25E2%2580%2593Scholes/TPPG5YWQ">Black–Scholes</a> option pricing model's <a href="/facts/Differential_equation/9MGbE3Ca">differential equation</a> can be transformed into the heat equation, and thus <a href="/facts/Numerical_methods/QT9cYa9z">numerical solutions</a> for <a href="/facts/Valuation_of_options/11Ft9oEA">option pricing</a> can be obtained with the Crank–Nicolson method. </p><p>The importance of this for finance is that option pricing problems, when extended beyond the standard assumptions (e.g. incorporating changing dividends), cannot be solved in closed form, but can be solved using this method. Note however, that 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 <a href="/facts/Vanilla_option/5gJFXV3X">vanilla options</a>, this results in oscillation in the <a href="/facts/Greeks_(finance)/v52IJIq5">gamma value</a> around the <a href="/facts/Strike_price/WzUb8XG6">strike price</a>. Therefore, special damping initialization steps are necessary (e.g., fully implicit finite difference method). </p> <h2 id="see-also">See also</h2> <ul><li><a href="/facts/Financial_mathematics/P9otEltQ">Financial mathematics</a></li> <li><a href="/facts/Trapezoidal_rule/KEg8cLGa">Trapezoidal rule</a></li></ul> <h2 id="external-links">External links</h2> <ul><li><a href="http://www3.nd.edu/~dbalsara/Numerical-PDE-Course/">Numerical PDE Techniques for Scientists and Engineers</a>, open access Lectures and Codes for Numerical PDEs</li> <li><a href="http://scicomp.stackexchange.com/questions/7399/how-to-discretize-the-advection-equation-using-the-crank-nicolson-methods">An example of how to apply and implement the Crank–Nicolson method for the Advection equation</a></li></ul> <h2 id="references">References</h2> <ol> <li id="fn:1"><p>Tuncer Cebeci (2002). Convective Heat Transfer. Springer. ISBN 0-9668461-4-1. <a href="0-9668461-4-1" target="_blank">0-9668461-4-1</a> <a href="#fnref:1" class="footnote-back-ref">↩</a></p></li> <li id="fn:2"><p>Crank, J.; Nicolson, P. (1947). "A practical method for numerical evaluation of solutions of partial differential equations of the heat conduction type". Mathematical Proceedings of the Cambridge Philosophical Society. 43 (1): 50–67. Bibcode:1947PCPS...43...50C. doi:10.1017/S0305004100023197. S2CID 16676040. <a href="/wiki/John_Crank" target="_blank">/wiki/John_Crank</a> <a href="#fnref:2" class="footnote-back-ref">↩</a></p></li> <li id="fn:3"><p>Thomas, J. W. (1995). Numerical Partial Differential Equations: Finite Difference Methods. Texts in Applied Mathematics. Vol. 22. Berlin, New York: Springer-Verlag. ISBN 978-0-387-97999-1.. Example 3.3.2 shows that Crank–Nicolson is unconditionally stable when applied to u t = a u x x {\displaystyle u_{t}=au_{xx}} . <a href="978-0-387-97999-1" target="_blank">978-0-387-97999-1</a> <a href="#fnref:3" class="footnote-back-ref">↩</a></p></li> <li id="fn:4"><p>"Multi-Dimensional Parabolic Problems" (PDF). Computer Science Department. RPI. Retrieved 29 May 2016. <a href="http://www.cs.rpi.edu/~flaherje/pdf/pde5.pdf" target="_blank">http://www.cs.rpi.edu/~flaherje/pdf/pde5.pdf</a> <a href="#fnref:4" class="footnote-back-ref">↩</a></p></li> <li id="fn:5"><p>Wilmott, P.; Howison, S.; Dewynne, J. (1995). The Mathematics of Financial Derivatives: A Student Introduction. Cambridge Univ. Press. ISBN 0-521-49789-2. The Mathematics of Financial Derivatives Wilmott. <a href="0-521-49789-2" target="_blank">0-521-49789-2</a> <a href="#fnref:5" class="footnote-back-ref">↩</a></p></li> </ol>