Crank–Nicolson method
In numerical analysis, the Crank–Nicolson method is a finite difference method used for numerically solving the heat equation and similar partial differential equations.
The method was developed by John Crank and Phyllis Nicolson in the 1940s.
[3] However, the approximate solutions can still contain (decaying) spurious oscillations if the ratio of time step
times the thermal diffusivity to the square of space step,
, is large (typically, larger than 1/2 per Von Neumann stability analysis).
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.
[citation needed] The Crank–Nicolson method is based on the trapezoidal rule, giving second-order convergence in time.
For linear equations, the trapezoidal rule is equivalent to the implicit midpoint method[citation needed]—the simplest example of a Gauss–Legendre implicit Runge–Kutta method—which also has the property of being a geometric integrator.
In many problems, especially linear diffusion, the algebraic problem is tridiagonal and may be efficiently solved with the tridiagonal matrix algorithm, which gives a fast
, Given that the terms on the right-hand side of the equation are known, this is a tridiagonal problem, so that
A quasilinear equation, such as (this is a minimalistic example and not general) 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
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.
chosen as constant), an advective component (which means that the system is evolving in space due to a velocity field), which we choose to be a constant
time) transforms each component of the PDE into the following: Now we create the following constants to simplify the algebra: and substitute (2), (3), (4), (5), (6), (7),
We express this as a linear system problem: where 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): 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: The d vector here is used to hold the boundary conditions.
In this example it is a 12×1 vector: To find the concentration at any time, one must iterate the following equation: 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 can be solved with the Crank–Nicolson discretization of assuming that a square grid is used, so that
This equation can be simplified somewhat by rearranging terms and using the CFL number 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 Solving such a linear system is costly.
Hence an alternating-direction implicit method 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 tridiagonal matrix algorithm to be solved.
One option is to use Newton's method to converge on the prediction, but this requires the computation of the Jacobian.
is the velocity of the system, then the Crank–Nicolson prediction will be a fixed point of the map
Even for high-dimensional systems, iteration of this map can converge surprisingly quickly.
[5] Particularly, the Black–Scholes option pricing model's differential equation can be transformed into the heat equation, and thus numerical solutions for option pricing can be obtained with the Crank–Nicolson method.
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 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).
