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}) – 2*h*_{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 2*N*_{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}) – 2*h*_{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.

Was this topic helpful?