Solves Poisson’s or Helmholtz’s equation on a two-dimensional rectangle using a fast Poisson solver based on the HODIE finite-difference scheme on a uniform mesh.

An array of size nx by ny containing the solution at the grid points.

Let c = coeff_u, ax = ax, bx = bx, ay = ay, by = by, nx = nx and ny = ny.

imsl_f_fast_poisson_2d is based on the code HFFT2D by Boisvert (1984). It solves the equation

on the rectangular domain (ax, bx) × (ay, by) with a user-specified combination of Dirichlet (solution prescribed), Neumann (first-derivative prescribed), or periodic boundary conditions. The sides are numbered clockwise, starting with the right side.

When c = 0 and only Neumann or periodic boundary conditions are prescribed, then any constant may be added to the solution to obtain another solution to the problem. In this case, the solution of minimum ∞-norm is returned.

The solution is computed using either a second-or fourth-order accurate finite-difference approximation of the continuous equation. The resulting system of linear algebraic equations is solved using fast Fourier transform techniques. The algorithm relies on the fact that nx - 1 is highly composite (the product of small primes). For details of the algorithm, see Boisvert (1984). If nx - 1 is highly composite then the execution time of imsl_f_fast_poisson_2d is proportional to nxnylog2nx. If evaluations of p(x, y) are inexpensive, then the difference in running time between order = 2 and order = 4 is small.

The grid spacing is the distance between the (uniformly spaced) grid lines. It is given by the formulas hx = (bx - ax)/(nx -1) and hy = (by - ay)/(ny - 1). The grid spacings in the x and y directions must be the same, i.e., nx and ny must be such that hx is equal to hy. Also, as noted above, nx and ny must be at least 4. To increase the speed of the fast Fourier transform, nx - 1 should be the product of small primes. Good choices are 17, 33, and 65.

If -coeff_u is nearly equal to an eigenvalue of the Laplacian with homogeneous boundary conditions, then the computed solution might have large errors.

On some platforms, imsl_f_fast_poisson_2d can evaluate the user-supplied function fcn in parallel. This is done only if the function imsl_omp_options is called to flag user-defined functions as thread-safe. A function is thread-safe if there are no dependencies between calls. Such dependencies are usually the result of writing to global or static variables.

In this example, the equation

with the boundary conditions

on the bottom side and

on the other three sides is solved. The domain is the rectangle [0, ¼] ×[0, ½]. The output of imsl_f_fast_poisson_2d is a 17 × 33 table of values. The functions imsl_f_spline_2d_value are used to print a different table of values.