Stability condition  forces one to use a rather small time step, which for integration over large times may take a considerable computational time. (With modern computers, this is rarely an issue when the PDE involves only one spatial coordinate, but for problems with two (or three) spatial coordinates, this issue is still significant.) Therefore, methods whose stability conditions would be less restrictive than  (or even absent) are of considerable interest. A famous example is the Crank-Nicolson (CN) method, obtained from the following discretization of :
This method is stable unconditionally (i.e., for any r), and its accuracy is O(r2) + O(h2), as we will explain below. Its stencil is shown in Figure 3 by the six gray circles. From  one can see that Um+ cannot be determined in isolation. Rather, one has to determine the vector of the unknowns on the entire (n + 1)th time level. Introducing notations
one can rewrite the CN method  as a matrix equation:
where I is the M x M unit matrix and r is defined in .
The number of arithmetic operations required to solve a generic M x M matrix equation (i.e., a system of M linear equations with M unknowns) is on the order of M3. Thus, if the spatial domain is discretized by, say, 100 mesh points, solving such a generic M x M matrix equation would take about one million operations per time step. This is many orders of magnitude greater than the operations count for method , which requires only about 4M operations per time step. However, the matrix (on the LHS) of the CN method  is not generic. Rather, it is tridiagonal, meaning that only its entries on the main diagonal and the two adjacent sub- and superdiagonals are nonzero (see ). For such matrices, it is possible to solve  using only about 8M operations, a tremendous improvement compared to the generic count of about M3. The corresponding solution method is called the Thomas algorithm (invented in the 1940s) and is discribed in most textbooks on numerical solution of PDEs or on numerical analysis. This method is also rather easy to program. Thus, the CN method offers a considerable computational time saving compared to the explicit method , because the time step in the CN method is no longer restricted by  to ensure stability.
However, the time step in the CN method still needs to be about, or smaller than, the mesh size h to ensure the accuracy 0(h2). This is especially important if the initial condition u0 has sharp corners or is discontinuous. In such a case, it has been found empirically that nonsmooth features of u0(x) 'diffuse away' (as they do in the exact solution) and do not contaminate the numerical solution when t < h/(3VD). Also, the smaller the r, the faster these nonsmooth features 'diffuse away' in the numerical solution. Alternatively, for a nonsmooth u0, one can use another method:
which has the same accuracy and the stability property as the CN method, but smoothes out nonsmooth initial conditions much more successfully (i.e., for larger r). (The derivation of this formula uses the same idea as the derivation of the two-step method  in Numerical Methods for Local Models. Method  is to be started by computing Ulm by, say, the CN method. To program the rest of , one moves all of the Un + 's on the LHS and all the other terms on the RHS and then rewrites this as a matrix equation analogous to . The matrix on the LHS is tridiagonal and hence can be inverted time-efficiently using the Thomas algorithm.
We now briefly explain why the accuracy of the CN method is 0(t2) + O(h2); this will show how  can be generalized for the more general equation . First, note that for any function u(t), different discretization schemes approximate the derivative u with different accuracy. For example, ujt + T ) — ujt ) u (t) =--+ 0(t)
as can be shown using the Taylor expansion of u near time t. That is, the approximation at point t is more accurate if one uses the data points located symmetrically on both sides of t. Note then that points t and t + t are located symmetrically about the point t + 1 /2t. Therefore, a formula similar to  approximates u at that point with accuracy 0(t2):
Therefore, the quotient appearing on the LHS of  approximates ii with accuracy 0(t2) about the virtual node at (xm, tn + (t/2)), shown in Figure 3 with a cross. Finally, similarly to  and  it can be shown that u(t -
Combining  with the second equation in , one has that the RHS of  approximates SXu at the same virtual node with accuracy O(h2) + 0(t2). Thus, the overall accuracy of  is 0(t2) (from the LHS) plus O(h2) + 0(t2) (from the RHS), resulting in the overall value 0(t2) + O(h2).
Was this article helpful?