Documentation |
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 N_{1} segments of length h_{1}, the second into N_{2} segments of length h_{2}, and (N_{1} + 1) by (N_{2} + 1) points are introduced on the regular grid thus defined. The triangles are all congruent with sides h_{1}, h_{2} 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 N_{2} – 1 diagonal blocks, here called T, are themselves tridiagonal (N_{1} – 1) by (N_{1} – 1) matrices, with 2(h_{1}/h_{2} + h_{2}/h_{1}) on the diagonal and –h_{2}/h_{1} on the subdiagonals. The subdiagonal blocks, here called I, are –h_{1}/h_{2} times the unit N_{1} – 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 (N_{1} – 1) by (N_{1} – 1) matrix with Sij= sin(πij/N_{1}). Then S^{–1}TS = Λ, where Λ is a diagonal matrix with diagonal entries 2(h_{1}/h_{2} + h_{2}/h_{1}) – 2h_{2}/h_{1} cos(πi/N_{1}). w = SΛ^{–1}S^{–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 2N_{1}.
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 N_{1} – 1 decoupled tridiagonal systems of size N_{2} – 1. Thus, a solution algorithm would look like
Divide F into N_{2} – 1 blocks of length N_{1} – 1, and perform an inverse discrete sine transform on each block.
Reorder the elements and solve N_{1} – 1 tridiagonal systems of size N_{2} – 1, with 2(h_{1}/h_{2} + h_{2}/h_{1}) – 2h_{2}/h_{1} cos(πi/N_{1}) on the diagonal, and –h_{1}/h_{2} on the subdiagonals.
Reverse the reordering, and perform N_{2} – 1 discrete sine transforms on the blocks of length N_{1} – 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 in poisolv. The discrete sine transform and the inverse discrete sine transform are computed by dst and idst, respectively.