d03edf
d03edf
© Numerical Algorithms Group, 2002.
Purpose
D03EDF Elliptic PDE, solution of finite difference equations by a
multigrid technique
Synopsis
[a,rhs,us,u,numit,ub,ifail] = d03edf(ngx,ngy,a,rhs<,ub,maxit,acc,iout,ifail>)
Description
D03EDF solves, by multigrid iteration, the seven-point scheme
6 7
A u +A u
i,j i-1,j+1 i,j i,j+1
3 4 5
+A u +A u +A u
i,j i-1,j i,j ij i,j i+1,j
1 2
+A u +A u =f ,
i,j i,j-1 i,j i+1,j-1 ij
i=1,2,...,n ; j=1,2,...,n ,
x y
which arises from the discretization of an elliptic partial
differential equation of the form
(alpha)(x,y)U +(beta)(x,y)U +(gamma)(x,y)U +(delta)(x,y)U
xx xy yy x
+(epsilon)(x,y)U +(phi)(x,y)U=(psi)(x,y)
y
and its boundary conditions, defined on a rectangular region.
This we write in matrix form as
Au=f
Systems of linear equations, matching the seven-point stencil
defined above, are solved by a multigrid iteration. An initial
estimate of the solution must be provided by the user. A zero
guess may be supplied if no better approximation is available.
A 'smoother' based on incomplete Crout decomposition is used to
eliminate the high frequency components of the error. A
restriction operator is then used to map the system on to a
sequence of coarser grids. The errors are then smoothed and
prolongated (mapped onto successively finer grids). When the
finest cycle is reached, the approximation to the solution is
corrected. The cycle is repeated for MAXIT iterations or until
the required accuracy, ACC, is reached.
D03EDF will automatically determine the number l of possible
coarse grids, 'levels' of the multigrid scheme, for a particular
problem. In other words, D03EDF determines the maximum integer l
so that n and n can be expressed in the form
x y
l-1 l-1
n =m2 +1, n =n2 +1, with m>=2 and n>=2.
x y
It should be noted that the rate of convergence improves
significantly with the number of levels used, so that n and n
x y
should be carefully chosen so that n -1 and n -1 have factors of
x y
1
the form 2 , with l as large as possible. For good convergence
the integer l should be at least 2.
D03EDF has been found to be robust in application, but being an
iterative method the problem of divergence can arise. For a
strictly diagonally dominant matrix A
4 -- k
|A |> > |A |
ij -- ij
k/=4
no such problem is foreseen. The diagonal dominance of A is not a
necessary condition, but should this condition be strongly
violated then divergence may occur. The quickest test is to try
the routine.
Parameters
d03edf
Required Input Arguments:
ngx integer
ngy integer
a (:,7) real
rhs (:) real
Optional Input Arguments: <Default>
ub (ngx*ngy) real zeros(ngx*ngy,1)
maxit integer 15
acc real sqrt(eps)
iout integer 0
ifail integer -1
Output Arguments:
a (:,7) real
rhs (:) real
us (:) real
u (:) real
numit integer
ub (ngx*ngy) real
ifail integer