While the general strategy of Partial Differential Equation Toolbox™ software is to use the MATLAB® built-in solvers for sparse systems, there are situations where faster solution algorithms are available. One such example is found when solving Poisson's equation
–Δu = f in Ω
with Dirichlet boundary conditions, where Ω is a rectangle.
For the fast solution algorithms to work, the mesh on the rectangle must be a regular mesh. In this context it means that the first side of the rectangle is divided into N1 segments of length h1, the second into N2 segments of length h2, and (N1 + 1) by (N2 + 1) points are introduced on the regular grid thus defined. The triangles are all congruent with sides h1, h2 and a right angle in between.
The Dirichlet boundary conditions are eliminated in the usual way, and the resulting problem for the interior nodes is Kv = F. If the interior nodes are numbered from left to right, and then from bottom to top, the K matrix is block tridiagonal. The N2 – 1 diagonal blocks, here called T, are themselves tridiagonal (N1 – 1) by (N1 – 1) matrices, with 2(h1/h2 + h2/h1) on the diagonal and –h2/h1 on the subdiagonals. The subdiagonal blocks, here called I, are –h1/h2 times the unit N1 – 1 matrix.
The key to the solution of the problem Kv = F is that the problem Tw = f is possible to solve using the discrete sine transform. Let S be the (N1 – 1) by (N1 – 1) matrix with Sij= sin(πij/N1). Then S–1TS = Λ, where Λ is a diagonal matrix with diagonal entries 2(h1/h2 + h2/h1) – 2h2/h1 cos(πi/N1). w = SΛ–1S–1 f, but multiplying with S is nothing more than taking the discrete sine transform, and multiplying with S–1 is the same as taking the inverse discrete sine transform. The discrete sine transform can be efficiently calculated using the fast Fourier transform on a sequence of length 2N1.
Solving Tw = f using the discrete sine transform would not be an advantage in itself, since the system is tridiagonal and should be solved as such. However, for the full system Ky = F, a transformation of the blocks in K turns it into N1 – 1 decoupled tridiagonal systems of size N2 – 1. Thus, a solution algorithm would look like
Divide F into N2 – 1 blocks of length N1 – 1, and perform an inverse discrete sine transform on each block.
Reorder the elements and solve N1 – 1 tridiagonal systems of size N2 – 1, with 2(h1/h2 + h2/h1) – 2h2/h1 cos(πi/N1) on the diagonal, and –h1/h2 on the subdiagonals.
Reverse the reordering, and perform N2 – 1 discrete sine transforms on the blocks of length N1 – 1.
When using a fast solver such as this one, time and memory are also saved since the matrix K in fact never has to be assembled. A drawback is that since the mesh has to be regular, it is impossible to do adaptive mesh refinement.
The fast elliptic solver for Poisson's equation is implemented
poisolv. The discrete sine transform and the
inverse discrete sine transform are computed by