View source for Crank-Nicolson method

One of the most popular methods for the numerical integration (cf. [[Integration, numerical|Integration, numerical]]) of diffusion problems, introduced by J. Crank and P. Nicolson [[#References|[a1]]] in 1947. They considered an implicit [[Finite difference|finite difference]] scheme to approximate the solution of a non-linear differential system of the type which arises in problems of heat flow.

In order to illustrate the main properties of the Crank–Nicolson method, consider the following initial-boundary value problem for the [[Heat equation|heat equation]]

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c1202601.png" /></td> </tr></table>

To approximate the solution <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c1202602.png" /> of this model problem on <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c1202603.png" />, one introduces the grid

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c1202604.png" /></td> </tr></table>

where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c1202605.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c1202606.png" /> are positive integers, and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c1202607.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c1202608.png" /> are the step sizes in space and time, respectively. One looks for approximations <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c1202609.png" /> to <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026010.png" />, and to this end one replaces derivatives by finite-difference formulas. Setting <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026011.png" />, then the Crank–Nicolson method for the considered problem takes the form

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026012.png" /></td> </tr></table>

<img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026013.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026014.png" />, with boundary conditions

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026015.png" /></td> </tr></table>

and numerical initial condition

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026016.png" /></td> </tr></table>

where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026017.png" /> is an approximation to <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026018.png" />. Note that each term in the difference equation can be interpreted as an approximation to the corresponding one in the differential equation at <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026019.png" />.

Whenever the theoretical solution <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026020.png" /> is sufficiently smooth, it can be proved, by [[Taylor series|Taylor series]] expansions that there exists a positive constant <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026021.png" />, which depends only on <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026022.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026023.png" />, such that the truncation error

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026024.png" /></td> </tr></table>

satisfies

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026025.png" /></td> </tr></table>

<img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026026.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026027.png" />. Thus, the Crank–Nicolson scheme is consistent of second order both in time and space.

The stability study ([[#References|[a2]]], [[#References|[a3]]]) in the discrete <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026028.png" />-norm can be made by Fourier analysis or by the energy method. One introduces the discrete operator

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026029.png" /></td> </tr></table>

<img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026030.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026031.png" />, where it is assumed that <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026032.png" />; and sets <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026033.png" />. There exists a positive constant <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026034.png" />, which is independent of <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026035.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026036.png" />, such that

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026037.png" /></td> </tr></table>

<img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026038.png" />. The stability estimate holds without any restriction on the step sizes <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026039.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026040.png" />; thus, the Crank–Nicolson method is said to be unconditionally stable. Convergence is derived by consistency and stability and, as <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026041.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026042.png" />, one finds

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026043.png" /></td> </tr></table>

<img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026044.png" />. Stability can also be established in the discrete <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026045.png" />-norm by means of an energy argument [[#References|[a3]]]. Denoting <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026046.png" />, and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026047.png" />, where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026048.png" />, then the following stability estimate holds:

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026049.png" /></td> </tr></table>

<img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026050.png" />. Therefore, if the theoretical solution <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026051.png" /> is sufficiently smooth, then

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026052.png" /></td> </tr></table>

<img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026053.png" />. Note that the above inequality implies a convergence estimate in the maximum norm.

An important question is to establish a maximum principle for the approximations obtained with the Crank–Nicolson method, similar to the one satisfied by the solutions of the heat equation. Related topics are monotonicity properties and, in particular, the non-negativity (or non-positivity) of the numerical approximations. Maximum-principle and monotonicity arguments can be used to derive convergence in the <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026054.png" />-norm. It can be proved ([[#References|[a2]]], [[#References|[a3]]]) that a choice of the step sizes such that <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026055.png" /> and the condition

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026056.png" /></td> </tr></table>

imply

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026057.png" /></td> </tr></table>

<img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026058.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026059.png" />. Note that if <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026060.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026061.png" />, then the following stability estimate in the maximum norm holds:

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026062.png" /></td> </tr></table>

where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026063.png" /> and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026064.png" />. This estimate is still valid with <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026065.png" /> if <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026066.png" /> (see [[#References|[a4]]]), and also it holds without any restriction between the step sizes but then a value <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026067.png" /> has to be accepted in the bound ([[#References|[a3]]], [[#References|[a5]]]).

From a computational point of view, the Crank–Nicolson method involves a tridiagonal linear system to be solved at each time step. This can be carried out efficiently by Gaussian elimination techniques [[#References|[a2]]]. Because of that and its accuracy and stability properties, the Crank–Nicolson method is a competitive algorithm for the numerical solution of one-dimensional problems for the heat equation.

The Crank–Nicolson method can be used for multi-dimensional problems as well. For example, in the integration of an homogeneous [[Dirichlet problem|Dirichlet problem]] in a rectangle for the heat equation, the scheme is still unconditionally stable and second-order accurate. Also, the system to be solved at each time step has a large and sparse matrix, but it does not have a tridiagonal form, so it is usually solved by iterative methods. The amount of work required to solve such a system is sufficiently large, so other numerical schemes are also taken into account here, such as alternating-direction implicit methods [[#References|[a6]]] or fractional-steps methods [[#References|[a7]]]. On the other hand, it should be noted that, for multi-dimensional problems in general domains, the finite-element method is better suited for the spatial discretization than the finite-difference method is.

The Crank–Nicolson method can be considered for the numerical solution of a wide variety of time-dependent partial differential equations. Consider the [[Abstract Cauchy problem|abstract Cauchy problem]]

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026068.png" /></td> </tr></table>

where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026069.png" /> represents a partial differential operator (cf. also [[Differential equation, partial|Differential equation, partial]]; [[Differential operator|Differential operator]]) which differentiates the unknown function <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026070.png" /> with respect to the space variable <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026071.png" /> in its space domain in <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026072.png" />, and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026073.png" /> may be a vector function. In the numerical integration of such problem, one can distinguish two stages: space discretization and time integration.

For the spatial discretization one can use finite differences, finite elements, spectral techniques, etc., and then a system of ordinary differential equations is obtained, which can be written as

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026074.png" /></td> </tr></table>

where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026075.png" /> is the space discretization parameter (the spatial grid size of a finite-difference or finite-element scheme, the inverse of the highest frequency of a spectral scheme, etc.) and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026076.png" /> is a suitable approximation to the theoretical initial condition <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026077.png" />.

The phrase  "Crank–Nicolson method"  is used to express that the time integration is carried out in a particular way. However, there is no agreement in the literature as to what time integrator is called the Crank–Nicolson method, and the phrase sometimes means the trapezoidal rule [[#References|[a8]]] or the implicit midpoint method [[#References|[a6]]]. Let <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026078.png" /> be the time step and introduce the discrete time levels <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026079.png" />, <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026080.png" />. If one uses the trapezoidal rule, then the full discrete method takes the form

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026081.png" /></td> </tr></table>

while when the implicit midpoint method is considered, one obtains

<table class="eq" style="width:100%;"> <tr><td valign="top" style="width:94%;text-align:center;"><img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026082.png" /></td> </tr></table>

where <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026083.png" />, and <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026084.png" /> is the approximation to <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026085.png" />. Both methods coincide for linear autonomous problems, but they are different in general. For instance, the midpoint rule is symplectic [[#References|[a9]]], while the trapezoidal rule does not satisfy this property.

The Crank–Nicolson method applied to several problems can be found in [[#References|[a8]]], [[#References|[a10]]], [[#References|[a11]]], [[#References|[a12]]], [[#References|[a13]]] and [[#References|[a14]]].

====References====
<table><TR><TD valign="top">[a1]</TD> <TD valign="top">  J. Crank,   P. Nicolson,   "A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type"  ''Proc. Cambridge Philos. Soc.'' , '''43'''  (1947)  pp. 50–67</TD></TR><TR><TD valign="top">[a2]</TD> <TD valign="top">  K.W. Morton,   D.F. Mayers,   "Numerical solution of partial differential equations" , Cambridge Univ. Press  (1994)</TD></TR><TR><TD valign="top">[a3]</TD> <TD valign="top">  V. Thomee,   "Finite difference methods for linear parabolic equations"  P.G. Ciarlet (ed.)  J.L. Lions (ed.) , ''Handbook of Numerical Analysis'' , '''1''' , Elsevier  (1990)  pp. 9–196</TD></TR><TR><TD valign="top">[a4]</TD> <TD valign="top">  J.F.B.M. Kraaijevanger,   "Maximum norm contractivity of discretization schemes for the heat equation"  ''Appl. Numer. Math.'' , '''9'''  (1992)  pp. 475–496</TD></TR><TR><TD valign="top">[a5]</TD> <TD valign="top">  C. Palencia,   "A stability result for sectorial operators in Banach spaces"  ''SIAM J. Numer. Anal.'' , '''30'''  (1993)  pp. 1373–1384</TD></TR><TR><TD valign="top">[a6]</TD> <TD valign="top">  A.R. Mitchell,   D.F. Griffiths,   "The finite difference method in partial differential equations" , Wiley  (1980)</TD></TR><TR><TD valign="top">[a7]</TD> <TD valign="top">  N.N. Yanenko,   "The method of fractional steps" , Springer  (1971)</TD></TR><TR><TD valign="top">[a8]</TD> <TD valign="top">  J.G. Verwer,   J.M. Sanz-Serna,   "Convergence of method of lines approximations to partial differential equations"  ''Computing'' , '''33'''  (1984)  pp. 297–313</TD></TR><TR><TD valign="top">[a9]</TD> <TD valign="top">  J.M. Sanz-Serna,   M.P. Calvo,   "Numerical Hamiltonian problems" , Chapman&amp;Hall  (1994)</TD></TR><TR><TD valign="top">[a10]</TD> <TD valign="top">  Y. Tourigny,   "Optimal <img align="absmiddle" border="0" src="https://www.encyclopediaofmath.org/legacyimages/c/c120/c120260/c12026086.png" /> estimates for two time-discrete Galerkin approximations of a nonlinear Schroedinger equation"  ''IMA J. Numer. Anal.'' , '''11'''  (1991)  pp. 509–523</TD></TR><TR><TD valign="top">[a11]</TD> <TD valign="top">  G.D. Akrivis,   "Finite difference discretization of the Kuramoto–Sivashinsky equation"  ''Numer. Math.'' , '''63'''  (1992)  pp. 1–11</TD></TR><TR><TD valign="top">[a12]</TD> <TD valign="top">  S.K. Chung,   S.N. Ha,   "Finite element Galerkin solutions for the Rosenau equation"  ''Appl. Anal.'' , '''54'''  (1994)  pp. 39–56</TD></TR><TR><TD valign="top">[a13]</TD> <TD valign="top">  I.M. Kuria,   P.E. Raad,   "An implicit multidomain spectral collocation method for stiff highly nonlinear fluid dynamics problems"  ''Comput. Meth. Appl. Mech. Eng.'' , '''120'''  (1995)  pp. 163–182</TD></TR><TR><TD valign="top">[a14]</TD> <TD valign="top">  G. Fairweather,   J.C. Lopez-Marcos,   "Galerkin methods for a semilinear parabolic problem with nonlocal boundary conditions"  ''Adv. Comput. Math.'' , '''6'''  (1996)  pp. 243–262</TD></TR></table>