FISHPACK
A Poisson Equation Solver

FISHPACK is a FORTRAN77 library which solves several forms of Poisson's equation, by John Adams, Paul Swarztrauber, Roland Sweet.

FISHPACK is a package of subroutines for solving separable partial differential equations in various coordinate systems. Such equations include the Laplace, Poisson and Helmholtz equations and have the form:


        Uxx + Uyy            = 0         (Laplace)
        Uxx + Uyy            = F(X,Y)    (Poisson)
        Uxx + Uyy + lambda*U = F(X,Y)    (Helmholtz)

in two dimensional cartesian coordinates.

FISHPACK is not limited to the 2D cartesian case. It can set up and solve the equations in coordinate systems including:

2D cartesian,
2D polar,
2D cylindrical,
3D spherical,
3D axisymmetric spherical (cylindrical),
3D cartesian.

The algorithm uses 5 point finite differences and an evenly spaced grid. Other routines are included which solve related problems in which the derivative terms have coefficient functions.

In the documentation, the words "standard" and "staggered" grid are used. The difference between these grids is simple. In the one dimensional case, a standard grid of size H on the interval (A,B) would be (A, A+H, A+2*H, ..., B-H, B) whereas a staggered grid would be (A+H/2, A+3H/2, ..., B-3H/2, B-H/2). Versions of 2D solvers are offered below using both kinds of grids. Depending on the boundary conditions, or other singularities near the boundary, one or the other type of grid may be preferred.

The version of FISHPACK presented here has been altered and adapted somewhat. To obtain a clean and correct copy of the original version, go to http://www.netlib.org/fishpack/index.html

Languages:

FISHPACK is available in a FORTRAN77 version.

Related Programs:

BETIS, a FORTRAN77 program which solves Laplace's equation in a 2D region using the boundary element method.

FEM_50, a MATLAB program which solves Laplace's equation in an arbitrary region using the finite element method.

FEM2D_POISSON_RECTANGLE, a FORTRAN77 program which solves Poisson's equation in a 2D rectangle using the finite element method.

FFTPACK5, a FORTRAN90 library which computes Fast Fourier Transforms. A version of this library is included in FISHPACK.

SERBA, a FORTRAN77 program which solves problems in planar elasticity using the boundary element method.

SLATEC, a FORTRAN90 library which includes a copy of FISHPACK.

Author:

John Adams, Paul Swarztrauber, Roland Sweet.

Reference:

Ulrich Schumann, Roland Sweet,
A direct method for the solution of Poisson's equation with Neumann boundary conditions on a staggered grid of arbitrary size,
Journal of Computational Physics,
Volume 20, 1976, pages 171-182.
Paul Swarztrauber,
A direct method for the discrete solution of separable elliptic equations,
SIAM Journal on Numerical Analysis,
Volume 11, 1974, pages 1136-1150.
Paul Swarztrauber, Roland Sweet,
Efficient FORTRAN Subprograms for the Solution of Elliptic Equations,
NCAR Technical Report TN/IA-109,
National Center for Atmospheric Research, 1975.
Roland Sweet,
A cyclic reduction algorithm for solving block tridiagonal systems of arbitrary dimensions,
SIAM Journal on Numerical Analysis,
Volume 14, September 1977, pages 706-720.

Source Code:

fishpack.f, the source code;
fishpack.sh, commands to compile the source code;

Examples and Tests:

FISHPACK_PRB runs many tests on the FISHPACK software.

fishpack_prb.f, a sample calling program;
fishpack_prb.sh, commands to compile, link and run the sample calling program;
fishpack_prb_output.txt, the output file.

List of Routines:

BCRH
BLKTR1 solves the linear system
BLKTRI solve linear system derived from a separable elliptic equation.
BSRH
CBLKTR is a complex version of BLKTRI.
CCMPB computes the roots of the b polynomials using routine
CHKPR4 checks the input parameters.
CHKPRM checks the input parameters for errors.
CHKSN4 checks if the PDE that SEPX4 must solve is a singular operator.
CHKSNG checks if the PDE that SEPELI must solve is a singular operator.
CMGNBN: complex generalized Buneman algorithm, linear equation solver.
CMPCSG computes required cosine values in ascending order.
CMPMRG merges two ascending strings of numbers.
CMPOSD solves Poisson's equation for Dirichlet boundary conditions.
CMPOSN solves Poisson's equation with Neumann boundary conditions.
CMPOSP solves poisson equation with periodic boundary conditions.
CMPTR3 solves a tridiagonal system.
CMPTRX solves a system of linear equations where the
COFX sets coefficients in the x-direction.
COFX4 sets coefficients in the x-direction.
COFY sets coefficients in y direction
COMPB computes the roots of the b polynomials using subroutine
COSGEN computes required cosine values in ascending order.
COSQB backward cosine quarter wave transform.
COSQB1 is a utility routine for COSQB.
COSQF forward cosine quarter wave transform.
COSQF1 is a utility routine for COSQF.
COSQI initializes the cosine quarter wave transform.
COST cosine transform.
COSTI initializes the cosine transform.
CPADD computes the eigenvalues of the periodic tridiagonal matrix
cproc applies a sequence of matrix operations to the vector x and
cprocp applies a sequence of matrix operations to the vector x and
cprod applies a sequence of matrix operations to the vector x and
cprodp applies a sequence of matrix operations to the vector x and
DEFE4 first approximates the truncation error given by
DEFER first approximates the truncation error given by
DX computes second order finite difference
DX4 computes second order finite difference
DY computes second order finite difference
DY4 computes second order finite difference
EPMACH computes an approximate machiine epsilon (accuracy)
FDUMP creates an error dump.
GENBUN: generalized Buneman algorithm, linear equation solver.
HSTCRT: solves the standard five-point finite difference approximation on a staggered grid to the Helmholtz equation in cartesian coordinates.
HSTCS1 is a utility routine for HSTCSP.
HSTCSP solves the standard five-point finite difference approximation on a staggered grid to the modified Helmholtz equation in spherical coordinates assuming axisymmetry (no dependence on longitude).
HSTCYL solves the standard five-point finite difference approximation on a staggered grid to the modified helmholtz equation in cylindrical coordinates. This two-dimensional modified Helmholtz equation results from the Fourier transform of a three-dimensional Poisson equation.
HSTPLR solves the standard five-point finite difference approximation on a staggered grid to the Helmholtz equation in polar coordinates.
HSTSSP solves the standard five-point finite difference approximation on a staggered grid to the Helmholtz equation in spherical coordinates and on the surface of the unit sphere (radius of 1).
HW3CRT solves the standard seven-point finite difference approximation to the Helmholtz equation in Cartesian coordinates.
HWSCRT solves the standard five-point finite difference approximation to the Helmholtz equation in Cartesian coordinates.
HWSCS1 is a utility routine for HWSCSP.
HWSCSP solves a finite difference approximation to the modified Helmholtz equation in spherical coordinates assuming axisymmetry (no dependence on longitude). This two dimensional modified Helmholtz equation results from the Fourier transform of the three dimensional Poisson equation.
HWSCYL solves a finite difference approximation to the Helmholtz equation in cylindrical coordinates. This modified Helmholtz equation results from the Fourier transform of the three-dimensional Poisson equation.
HWSPLR solves a finite difference approximation to the Helmholtz equation in polar coordinates.
HWSSS1 is a utility routine for HWSSSP.
HWSSSP solves a finite difference approximation to the Helmholtz equation in spherical coordinates and on the surface of the unit sphere (radius of 1).
INDXB indexes the first root of the B(I,IR) polynomial.
INDXC
INXCA
INXCB
INXCC
J4SAVE sets or gets variables needed by the error handler.
MERGE merges two ascending strings of numbers in the array TCOS.
MINSO4 orthogonalizes the array usol with respect to
MINSOL orthogonalizes the array usol with respect to
ORTHO4 orthogonalizes the array usol with respect to
ORTHOG orthogonalizes the array usol with respect to
PGSF
PIMACH supplies the value of the constant pi correct to
POIS3D solves a special set of linear equations.
POISD2 solves Poisson's equation for Dirichlet boundary conditions.
POISN2 solves Poisson's equation with Neumann boundary conditions.
POISP2 solves Poisson's equation with periodic boundary conditions.
POISTG solves a special set of linear equations.
POS3D1
POSTG2 solves Poisson's equation on a staggered grid.
PPADD computes the eigenvalues of the periodic tridiagonal matrix
PROC applies a sequence of matrix operations to the vector x and
PROCP applies a sequence of matrix operations to the vector x and
PROD applies a sequence of matrix operations to the vector x and
PRODP applies a sequence of matrix operations to the vector x
RADB2 - backward Fourier transform, radix 2.
RADB3 - backward Fourier transform, radix 3.
RADB4 - backward Fourier transform, radix 4.
RADB5 - backward Fourier transform, radix 5.
RADBG - backward Fourier transform, general radix.
RADF2 - forward Fourier transform, radix 2.
RADF3 - forward Fourier transform, radix 3.
RADF4 - forward Fourier transform, radix 4.
RADF5 - forward Fourier transform, radix 5.
RADFG - forward Fourier transform, general radix.
RFFTB - backward Fourier transform.
RFFTB1
RFFTF - forward Fourier transform.
RFFTF1
RFFTI - initialized Fourier transform.
RFFTI1
SEPELI 2D general separable elliptic problem, second or fourth order scheme.
SEPX4 2D restricted separable elliptic problem, second or fourth order scheme.
SINQB backward sine quarter wave transform.
SINQF forward sine quarter wave transform.
SINQI initializes the sine quarter wave transform.
SINT: the sine transform.
SINTI initializes the sine transform.
SPELI4 sets up vectors and arrays for input to BLKTRI
SPELIP sets up vectors and arrays for input to BLKTRI
STORE forces its argument to be stored.
TEVLC
TEVLS finds the eigenvalues of a symmetric tridiagonal matrix.
TRI3
TRID
TRIS4 solves for a non-zero eigenvector corresponding
TRISP solves for a non-zero eigenvector corresponding
TRIX solves a system of linear equations where the
XERCNT allows the user to control error handling.
XERHLT aborts the program and prints an error message.
XERMAX sets the maximum number of appearances of an error message.
XERMSG processes an error message.
XERPRN prints an error message.
XERSVE records that an error has occurred.
XGETUA returns error unit numbers.
XSETF sets the error control flag.
XSETUN sets the error message output unit.

You can go up one level to the FORTRAN77 source codes.

Last revised on 15 October 2012.

FISHPACK A Poisson Equation Solver