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.
Synopsis
#include <imsl.h>
float *imsl_f_fast_poisson_2d (float rhs_pde(), float rhs_bc(), float coeff_u, int nx, int ny, float ax, float bx, float ay, float by, Imsl_bc_type bc_type[], ..., 0)
The type double function is imsl_d_fast_poisson_2d.
Required Arguments
float rhs_pde
(float
x,
float
y)
User-supplied function to evaluate the right-hand side of the
partial differential equation at x and y.
float
rhs_bc(Imsl_pde_side side, float x, float y)
User-supplied
function to evaluate the right-hand side of the boundary conditions, on side
side, at x and y. The value of side will be one of
the following: IMSL_RIGHT, IMSL_BOTTOM, IMSL_LEFT, or IMSL_TOP.
float coeff_u
(Input)
Value of the coefficient of u in the differential
equation.
int nx
(Input)
Number of grid lines in the x-direction. nx must be at least 4.
See the “Description” section for further
restrictions on nx.
int ny
(Input)
Number of grid lines in the y-direction. ny must be at least 4.
See the “Description” section for further
restrictions on ny.
float ax
(Input)
The value of x along the left side of the domain.
float bx
(Input)
The value of x along the right side of the domain.
float ay
(Input)
The value of y along the bottom of the domain.
float by
(Input)
The value of y along the top of the domain.
Imsl_bc_type
bc_type[4] (Input)
Array of size 4 indicating the type of
boundary condition on each side of the domain or that the solution is periodic.
The sides are numbered as follows:
Side
Location
IMSL_RIGHT_SIDE(0)
x = bx
IMSL_BOTTOM_SIDE(1)
y = ay
IMSL_LEFT_SIDE(2)
x = ax
IMSL_TOP_SIDE(3)
y = by
The three possible boundary condition types are as follows:
Type
Condition
IMSL_DIRICHLET
Value of u is given.
IMSL_NEUMANN
Value of du/dx is given (on the right or
left sides) or du/dy (on the bottom or
top of the domain).
IMSL_PERIODIC
Periodic.
Synopsis with Optional Arguments
#include <imsl.h>
float
*imsl_f_fast_poisson_2d (float rhs_pde(), float rhs_bc(),
float coeff_u, int nx, int ny, float ax,
float bx, float ay, float by, Imsl_bc_type
bc_type[],
IMSL_RETURN_USER, float u_user[],
IMSL_ORDER, int
order,
IMSL_RHS_PDE_W_DATA, float rsh_pde (), void
*data,
IMSL_RHS_BC_W_DATA, float rsh_bc (), void
*data,
0)
Optional Arguments
IMSL_RETURN_USER, float u_user[]
(Output)
User-supplied array of size nx by ny containing solution
at the grid points.
IMSL_ORDER, int order
(Input)
Order of accuracy of the finite-difference approximation. It can be
either 2 or 4.
Default: order = 4
IMSL_RSH_PDE_W_DATA, float
rhs_pde (float x,
float y,
void *data), void
*data,
(Input)
User-supplied function to evaluate the right-hand side of the partial
differential equation at x and y, which also accepts
a pointer to data that is supplied by the user. data is a pointer to the
data to be passed to the user-supplied function. See the Introduction, Passing Data to
User-Supplied Functions at the beginning of this manual for more details.
IMSL_RSH_BC_W_DATA,
float
rhs_bc(Imsl_pde_side side,
float x,
float y,
void *data) , void
*data,
(Input)
User-supplied function to evaluate right-hand side of the boundary
conditions, which also accepts a pointer to data that is supplied by the
user. data is a pointer to the data to be passed to the user-supplied
function. See the Introduction, Passing Data to
User-Supplied Functions at the beginning of this manual for more details.
Description
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 nxny log2 nx. 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.
Example
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.
#include <imsl.h>
#include <math.h>
main()
{
float rhs_pde(float, float);
float rhs_bc(Imsl_pde_side, float, float);
int nx =
17;
int nxtabl = 5;
int ny = 33;
int nytabl = 5;
int
i;
int j;
Imsl_f_spline *sp;
Imsl_bc_type bc_type[4];
float ax, ay, bx,
by;
float x, y, xdata[17], ydata[33];
float coefu, *u;
float u_table;
float abs_error;
/* Set rectangle size */
ax = 0.0;
bx = 0.25;
ay = 0.0;
by = 0.50;
/* Set boundary conditions */
bc_type[IMSL_RIGHT_SIDE] = IMSL_DIRICHLET_BC;
bc_type[IMSL_BOTTOM_SIDE] = IMSL_NEUMANN_BC;
bc_type[IMSL_LEFT_SIDE] = IMSL_DIRICHLET_BC;
bc_type[IMSL_TOP_SIDE] = IMSL_DIRICHLET_BC;
/* Coefficient of u */
coefu = 3.0;
/* Solve the PDE */
u =
imsl_f_fast_poisson_2d(rhs_pde, rhs_bc, coefu, nx, ny,
ax, bx, ay, by, bc_type, 0);
/* Set up for interpolation */
for (i = 0; i
< nx; i++)
xdata[i] = ax + (bx - ax) * (float) i / (float) (nx - 1);
for (i = 0; i
< ny; i++)
ydata[i] = ay + (by - ay) * (float) i / (float) (ny - 1);
/* Compute interpolant */
sp =
imsl_f_spline_2d_interp(nx, xdata, ny, ydata, u, 0);
printf("
x
y
u error\n\n");
for (i = 0; i < nxtabl; i++)
for (j = 0; j < nytabl; j++) {
x = ax + (bx - ax) * (float) j / (float) (nxtabl -
1);
y = ay + (by - ay) * (float) i / (float) (nytabl -
1);
u_table = imsl_f_spline_2d_value(x, y, sp, 0);
abs_error = fabs(u_table - sin(x + 2.0 * y) -
exp(2.0 * x + 3.0 * y));
/* Print computed answer and absolute on
nxtabl by nytabl grid */
printf(" %6.4f %6.4f
%6.4f %8.2e\n",
x, y, u_table, abs_error);
}
}
float rhs_pde(float x, float y)
{
/* Define the right side of the PDE */
return (-2.0 *
sin(x + 2.0 * y) + 16.0 * exp(2.0 * x + 3.0 * y));
}
float rhs_bc(Imsl_pde_side side, float x, float y)
{
/* Define the boundary conditions */
if (side ==
IMSL_BOTTOM_SIDE)
return (2.0 * cos(x + 2.0 * y) + 3.0 * exp(2.0 * x + 3.0 *
y));
else
return (sin(x + 2.0 * y) + exp(2.0 * x + 3.0 * y));
}
Output
x y u error
0.0000 0.0000 1.0000 0.00e+00
0.0625 0.0000 1.1956 5.12e-06
0.1250 0.0000 1.4087 7.19e-06
0.1875 0.0000 1.6414 5.10e-06
0.2500 0.0000 1.8961 8.67e-08
0.0000 0.1250 1.7024 1.73e-07
0.0625 0.1250 1.9562 6.39e-06
0.1250 0.1250 2.2345 9.50e-06
0.1875 0.1250 2.5407 6.36e-06
0.2500 0.1250 2.8783 1.66e-07
0.0000 0.2500 2.5964 2.60e-07
0.0625 0.2500 2.9322 9.25e-06
0.1250 0.2500 3.3034 1.34e-05
0.1875 0.2500 3.7148 9.27e-06
0.2500 0.2500 4.1720 9.40e-08
0.0000 0.3750 3.7619 4.84e-07
0.0625 0.3750 4.2163 9.16e-06
0.1250 0.3750 4.7226 1.36e-05
0.1875 0.3750 5.2878 9.44e-06
0.2500 0.3750 5.9199 5.72e-07
0.0000 0.5000 5.3232 5.93e-07
0.0625 0.5000 5.9520 9.84e-07
0.1250 0.5000 6.6569 1.34e-06
0.1875 0.5000 7.4483 4.55e-07
|
Visual Numerics, Inc. PHONE: 713.784.3131 FAX:713.781.9260 |