Such problems are considerably more difficult than those with one dimension. First, if the spatial region does not have rectangular or circular symmetry, the numerical modeling of the problem requires the use of finite-element methods. In circular regions, the problem needs to be reformulated in polar coordinates. Below we consider
Figure 5 The stencil of method .
Figure 5 The stencil of method .
only methods for rectangular regions. As earlier, we explain the main issues using the two-dimensional Heat equation
with conditions along the boundary of the rectangle 0 < x < X, 0 < y < Y which will be specified in each particular case. Generalizations to more realistic equations, such as a two-dimensional counterpart of , can be made along the lines considered above. Without loss of generality we assume that the mesh size along the x- and y-directions is the same and equals h; the numerical solution at node (xm = mh, yi = ih, tn = nr) is denoted U"mi, 0 < m < M, 0 < l < L. Two time levels for such a discretization are shown in Figure 4, with the boundary being shown by open circles.
First, the explicit Euler method for  is
where r is defined in  with D = 1. The stencil of this method is shown in Figure 5. The accuracy of this method is O(r) + O(h2) and its stability condition is r < 1/4
Any of the Dirichlet, Neumann, or mixed boundary conditions can be used with this method similarly to how it was done in one dimension. In addition, if either periodic or zero-flux (a particular case of Neumann) boundary conditions can be used, one can use two-dimensional built-in commands to approximate the RHS of ; in Matlab such commands are fft2 and dct2. The stability condition in this case is more restrictive:
A naive generalization of the CN method has the form n+1 ,l
To represent this as systems of linear equations, one arranges the (L — 1) x (M — 1) two-dimensional array Umj into a (L — 1) (M — 1)-component vector using the lexicographical order, an example of which is shown in the lower time level of Figure 4. The reason we called this generalization 'naive' is because although it has the same accuracy and stability as the CN method in one dimension, the coefficient matrix on the LHS of  is not tridiagonal. Therefore, it cannot be inverted time-efficiently. There are two ways out of this difficulty. One and the standard in the community of computational scientists is to use any of alternating direction implicit (ADI) methods; two well-known representatives of this family are the Peaceman-Rachford and DouglasRachford methods. Roughly speaking, these are special cases of operator-splitting methods, which perform the operations dx and 8j, in sequence. These methods are optimally time-efficient. However, both their programming and treatment of boundary conditions are nontrivial and are better be left to experts. The other way is to attempt to invert the nontridiagonal matrix in question using the fact that it is sparse (i.e., has nonzero entries only on five (not adjacent) diagonals). Matlab has special algorithms for dealing with such matrices (see the help entries for sparse and spdiags) which speed up the computations compared to the case of nonsparse matrices. The time efficiency of such an approach is suboptimal, but it may be a worthwhile sacrifice for the ease of the programming. Alternatives to these include using fft2 and dct2 when boundary conditions allow, or returning to the explicit method, also sacrificing the computational speed for the simplicity of coding.
Was this article helpful?