Note: This page has been translated by MathWorks. Please click here

To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

To view all translated materals including this page, select Japan from the country navigator on the bottom of this page.

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?