pde_1d_mg
Solves a system of one-dimensional time-dependent partial differential equations using a moving-grid interface.
Synopsis
#include <imsl.h>
void imsl_f_pde_1d_mg_mgr (int task, void **state,…, 0)
void imsl_f_pde_1d_mg (int npdes, int ngrids, float *t, float tend float u[], float xl, float xr, void *state, void pde_systems(), void boundary_conditions(), …, 0)
The void functions imsl_d_pde_1d_mg_mgr and imsl_d_pde_1d_mg are for double type arithmetic accuracy.
The function imsl_f_pde_1d_mg_mgr is used to initialize and reset the problem, and the function imsl_f_pde_1d_mg is the integrator. The descriptions of both functions are provided below.
|
NOTE: The integrator is provided with single or double precision arithmetic. We recommend using the double precision interface imsl_d_pde_1d_mg. |
Required Arguments for imsl_f_pde_1d_mg_mgr
int task (Input)
This function must be called with task set to IMSL_PDE_INITIALIZE to set up for solving a system and with task equal to IMSL_PDE_RESET to clean up after it has been solved. These values for task are defined in the include file, imsl.h.
void **state (Input/Output)
The current state of the PDE solution is held in a structure pointed to by state. It cannot be directly manipulated.
Required Arguments for imsl_f_pde_1d_mg
int npdes (Input)
The number of differential equations.
int ngrids (Input)
The number of spatial grid/mesh points, including the boundary points 
and 
.
float *t (Input/Output)
On input, t is the initial independent variable value. On output, t is replaced by tend, unless error conditions arise. This is first set to the value of the independent variable 
where the integration of 
begins. It is set to the value tend on return.
float tend (Input)
Mathematical value of 
where the integration of 
ends. Note: Starting values of t < tend imply integration in the forward direction, while values of t > tend imply integration in the backward direction. Either direction is permitted.
float u[] (Input/Output)
Array of size npdes+1 by ngrids. On input, the first npdes rows contain initial values for all components of the system at the equally spaced grid of values. It is not required to define the grid values in the last row of u. On output u[] contains the approximate solution value U i(xj(tend), tend) at array location u[i×ngrids+j]. The grid value x j(tend) is in location u[(npdes*ngrids) +j]. Normally the grid values are equally spaced as the integration starts. Variable grid values can be provided by defining them as output from the user function initial_conditions supplied by either imsl_f_pde_1d_mg_mgr’s IMSL_INITIAL_CONDITIONS, or IMSL_INITIAL_CONDITIONS_W_DATA optional arguments.
float xl (Input)
Lower grid boundary, 
.
float xr (Input)
Upper grid boundary, 
.
void *state (Input/Output)
The current state of the solution is held in a structure pointed to by state. It must be initialized by a call to imsl_f_pde_1d_mg_mgr. It cannot be directly manipulated.
void pde_systems(float t, float x, int npdes, int ngrids, float *full_u, float *grid_u, float *dudx, float *c, float *q, float *r, int *ires) (Input)
A user-supplied function to evaluate the differential equation, as expressed in Equation 2. Each application requires a function specifically designed for the task, and this function is normally written by the user of the integrator.
Evaluate the terms of the system of Equation 2. A default value of 
is assumed, but this can be changed to one of the choices, 
. Use the optional arguments IMSL_CART_COORDINATES, IMSL_CYL_COORDINATES, IMSL_SPH_COORDIANTES for the respective values 
. Return the values in the arrays as indicated:
If any of the functions cannot be evaluated, set ires=3. Otherwise, do not change the value of ires.
void boundary_conditions (float t, float *beta, float *gamma, float *full_u, float *grid_u, float *dudx, int npdes, int grids, int left, int *ires) (Input)
User-supplied function to supply the boundary conditions, as expressed in Equation 2.
The value 
, and the flag left=1 for 
. The flag has the value left=0 for 
. If any of the functions cannot be evaluated, set ires=3. Otherwise, do not change the value of ires.
Synopsis with Optional Arguments for imsl_f_pde_1d_mg_mgr
#include <imsl.h>
void imsl_f_pde_1d_mg_mgr (int task, void **state,
IMSL_CART_COORDINATES, or
IMSL_CYL_COORDINATES, or
IMSL_SPH_COORDINATES,
IMSL_TIME_SMOOTHING, float tau,
IMSL_SPATIAL_SMOOTHING, float kappa,
IMSL_MONITOR_REGULARIZING, float alpha,
IMSL_MAX_BDF_ORDER, int max_bdf_order,
IMSL_USER_FACTOR_SOLVE, int fac(), void sol(),
IMSL_USER_FACTOR_SOLVE_W_DATA, int fac(), void sol(), void data,
IMSL_INITIAL_CONDITIONS, void initial_conditions()
IMSL_INITIAL_CONDITIONS_W_DATA, void initial_conditions(), void data,
0)
Optional Arguments
IMSL_CART_COORDINATES, or
IMSL_CYL_COORDINATES, or
IMSL_SPH_COORDINATES
IMSL_CART_COORDINATES specifies cartesian coordinates, where 
in Equation 2. IMSL_CYL_COORDINATES specifies cylindrical or polar coordinates, where 
in Equation 2. IMSL_SPH_COORDINATES specifies spherical coordinates, where 
in Equation 2.
Default: IMSL_CART_COORDINATES
IMSL_TIME_SMOOTHING, float tau, (Input)
Resets the value of the parameter 
, described above.
Default: 
.
IMSL_SPATIAL_SMOOTHING, float kappa, (Input)
Resets the value of the parameter 
, described above.
Default: 
.
IMSL_MONITOR_REGULARIZING, float alpha, (Input)
Resets the value of the parameter 
, described above.
Default: 
.
IMSL_MAX_BDF_ORDER, int max_bdf_order, (Input)
Resets the maximum order for the bdf formulas used in imsl_f_dea_petzold_gear. The new value can be any integer between 1 and 5. Some problems benefit by making this change. The default value of max_bdf_order was chosen because imsl_f_dea_petzold_gear may cycle on its selection of order and step-size with max_bdf_order higher than value 2.
Default: max_bdf_order=2.
IMSL_USER_FACTOR_SOLVE, int fac(int neq, int iband, float *a), void sol(int neq, int iband, float *g, float *y) (Input)
User-supplied functions to factor A, and solve the system AΔy = Δg. Use of this optional argument allows for handling the factorization and solution steps in a problem‑specific manner. If successful fac should return 0, if unsuccessful, fac should return a non-zero value. See Example 5 - A Flame Propagation Model for sample usage of this optional argument.
IMSL_USER_FACTOR_SOLVE_W_DATA, int fac(int neq,int iband,float *a, void *data), void sol(int neq,int iband,float *g, float *y, void *data), void *data (Input)
User-supplied functions to factor A, and solve the system AΔy = Δg.The argument data is a pointer to the data that is passed to the user-supplied function.
IMSL_INITIAL_CONDITIONS, void initial_conditions(int npdes, int ngrids, float *u) (Input)
User-supplied function to supply the initial values for the system at the starting independent variable value t. This function can also provide a non-uniform grid at the initial value. Here npdes is the number of differential equations, ngrids is the number of grid points, and u is an array of size npdes+1 by ngrids, containing the approximate solution value 
in location u[i×ngrids+j]. The grid values are equally spaced on input, but can be updated as desired, provided the values are increasing. Update the grid values in array locations u[(npdes × ngrids) +j], where 0 ≤ j ≤ ngrids.
IMSL_INITIAL_CONDITIONS_W_DATA, void initial_conditions(int npdes, int ngrids, float *u, float *grid, void *data), void *data (Input)
User-supplied function to supply the initial values for the system at the starting independent variable value t. This function can also provide a non-uniform grid at the initial value. The argument data is a pointer to the data that is passed to the user-supplied function.
Synopsis with Optional Arguments for imsl_f_pde_1d_mg
#include <imsl.h>
void imsl_f_pde_1d_mg (int npdes ,int ngrids ,float *t ,float tend,float u[] ,float xl ,float xr ,void *state ,void pde_systems() ,void boundary_conditions(),
IMSL_RELATIVE_TOLERANCE ,float rtol,
IMSL_ABSOLUTE_TOLERANCE ,float atol,
IMSL_PDE_SYS_W_DATA ,void pde_systems() ,void *data,
IMSL_BOUNDARY_COND_W_DATA ,void boundary_conditions() ,void *data,
0)
Optional Arguments
IMSL_RELATIVE_TOLERANCE, float rtol, (Input)
This option resets the value of the relative accuracy parameter used in imsl_f_dea_petzold_gear.
Default: rtol=1.0E-2 for single precision, rtol=1.0E-4 for double precision.
IMSL_ABSOLUTE_TOLERANCE, float atol, (Input)
This option resets the value of the absolute accuracy parameter used in imsl_f_dea_petzold_gear.
Default: atol=1E-2 for single precision, atol=1E-4 for double precision.
IMSL_PDE_SYS_W_DATA, void pde_systems(float t, float x, int npdes, int ngrids, float *full_u, float *grid_u, float *dudx, float *c, float *q, float *r, int *ires, void *data),void *data (Input)
User-supplied function to evaluate the differential equation, as expressed in Equation 2. The argument data is a pointer to the data that is passed to the user-supplied function.
IMSL_BOUNDARY_COND_W_DATA, void boundary_conditions (float t,float *beta, float *gamma, float *full_u, float *grid_u, float *dudx, int npdes, int ngrids, int left,int *ires, void *data), void *data (Input)
User-supplied function to supply the boundary conditions, as expressed in Equation 2. The argument data is a pointer to the data that is passed to the user-supplied function.
Examples
Due to its importance and the complexity of its interface, function imsl_f_pde_1d_mg is presented with several examples. Many of the program features are exercised. The problems complete without any change to the optional arguments, except where these changes are required to describe or to solve the problem.
In many applications the solution to a PDE is used as an auxiliary variable, perhaps as part of a larger design or simulation process. The truncation error of the approximate solution is commensurate with piece-wise linear interpolation on the grid of values, at each output point. To show that the solution is reasonable, a graphical display is revealing and helpful. We have not provided graphical output as part of our documentation, but users may already have the Rogue Wave, Inc. product, PV‑WAVE, which is not included with IMSL C Numerical Library. Examples 1 through 8 write results in files pde_ex0#.out that can be visualized with PV‑WAVE. We provide a script of commands, pde_1d_mg_plot.pro, for viewing the solutions. This is listed below. The grid of values and each consecutive solution component is displayed in separate plotting windows. The script and data files written by examples 1-8 on a SUN-SPARC system are in the directory for IMSL C Numerical Library examples. When executing PV_WAVE, use the command line
pde_1d_mg_plot,filename=’pde_ex0#.out’
to view the output of a particular example. The symbol ‘#’ will be one of the choices 1,2,…,8. However, it is not necessary to have PV_WAVE installed to execute the examples.
To view the code, see Code for PV-WAVE Plotting.
Example 1 - Electrodynamics Model
This example is from Blom and Zegeling (1994). The system is
We make the connection between the model problem statement and the example:
The boundary conditions are
This is a non-linear problem with sharply changing conditions near 
. The default settings of integration parameters allow the problem to be solved. The use of imsl_f_pde_1d_mg requires two subroutines provided by the user to describe the differential equations, and boundary conditions.
#include <stdio.h>
#include <math.h>
#include <imsl.h>
/* prototypes */
static void initial_conditions (int npdes, int ngrids, double u[]);
static void pde_systems (double t, double x, int npdes, int ngrids,
double full_u[], double grid_u[], double dudx[], double *c,
double q[], double r[], int *ires);
static void boundary_conditions (double t, double beta[], double gamma[],
double full_u[], double grid_u[], double dudx[], int npdes,
int ngrids, int left, int *ires);
#define MIN(X,Y) (X<Y)?X:Y
#define NPDE 2
#define NFRAMES 5
#define N 51
#define U(I_,J_) u[I_ * ngrids + J_]
int main ()
{
char *state = NULL;
int i, j;
double u[(NPDE + 1) * N];
double t0 = 0.0, tout;
double delta_t = 10.0, tend = 4.0;
int npdes = NPDE, ngrids = N;
double xl = 0.0, xr = 1.0;
FILE *file1;
file1 = fopen ("pde_ex01.out", "w");
imsl_output_file (IMSL_SET_OUTPUT_FILE, file1, 0);
fprintf (file1, " %d\t%d\t%d", npdes, ngrids, NFRAMES);
fprintf (file1, "\t%f\t%f\t%f\t%f\n", xl, xr, t0, tend);
/* initialize u */
initial_conditions (npdes, ngrids, u);
imsl_d_pde_1d_mg_mgr (IMSL_PDE_INITIALIZE, &state, 0);
tout = 1e-3;
do
{
imsl_d_pde_1d_mg (npdes, ngrids, &t0, tout, u, xl, xr, state,
pde_systems, boundary_conditions, 0);
fprintf (file1, "%f\n", tout);
for (i = 0; i < npdes + 1; i++)
{
for (j = 0; j < ngrids; j++)
{
fprintf (file1, "%16.10f ", U (i, j));
if (((j + 1) % 4) == 0)
fprintf (file1, "\n");
}
fprintf (file1, "\n");
}
t0 = tout;
tout = tout * delta_t;
tout = MIN (tout, tend);
}
while (t0 < tend);
imsl_d_pde_1d_mg_mgr (IMSL_PDE_RESET, &state, 0);
#undef MIN
#undef NPDE
#undef NFRAMES
#undef N
#undef U
}
static void
initial_conditions (int npdes, int ngrids, double u[])
{
#define U(I_,J_) u[I_ * ngrids + J_]
int i;
for (i = 0; i < ngrids; i++)
{
U (0, i) = 1.0;
U (1, i) = 0.0;
}
#undef U
}
static void
pde_systems (double t, double x, int npdes, int ngrids,
double full_u[], double grid_u[], double dudx[], double *c,
double q[], double r[], int *ires)
{
#define C(I_,J_) c[I_ * npdes + J_]
double z;
static double eps = 0.143;
static double eta = 17.19;
static double pp = 0.1743;
C (0, 0) = 1.0;
C (0, 1) = 0.0;
C (1, 0) = 0.0;
C (1, 1) = 1.0;
r[0] = pp * dudx[0] * eps;
r[1] = pp * dudx[1];
z = eta * (grid_u[0] - grid_u[1]) / 3.0;
q[0] = exp (z) - exp (-2.0 * z);
q[1] = -q[0];
return;
#undef C
}
static void
boundary_conditions (double t, double beta[], double gamma[],
double full_u[], double grid_u[], double dudx[],
int ngrids, int npdes, int left, int *ires)
{
if (left)
{
beta[0] = 1.0;
beta[1] = 0.0;
gamma[0] = 0.0;
gamma[1] = grid_u[1];
}
else
{
beta[0] = 0.0;
beta[1] = 1.0;
gamma[0] = grid_u[0] - 1.0;
gamma[1] = 0.0;
}
return;
}
Example 2 - Inviscid Flow on a Plate
This example is a first order system from Pennington and Berzins, (1994). The equations are
Following elimination of w, there remain 
differential equations. The variable 
is not time, but a second space variable. The integration goes from 
to 
. It is necessary to truncate the variable 
at a finite value, say 
. In terms of the integrator, the system is defined by letting 
and
The boundary conditions are satisfied by
We use 
grid points and output the solution at steps of 
.
This is a non-linear boundary layer problem with sharply changing conditions near 
. The problem statement was modified so that boundary conditions are continuous near 
. Without this change the underlying integration software, imsl_f_dea_petzold_gear, cannot solve the problem. The continuous blending function 
is arbitrary and artfully chosen. This is a mathematical change to the problem, required because of the stated discontinuity at 
. Reverse communication is used for the problem data. No additional user-written subroutines are required when using reverse communication. We also have chosen 10 of the initial grid points to be concentrated near 
, anticipating rapid change in the solution near that point. Optional changes are made to use a pure absolute error tolerance and non-zero time-smoothing.
#include <stdio.h>
#include <math.h>
#include <imsl.h>
/* prototypes */
static void initial_conditions (int npdes, int ngrids, double u[]);
static void pde_systems (double t, double x, int npdes, int ngrids,
double full_u[], double grid_u[],
double dudx[], double *c, double q[],
double r[], int *ires);
static void boundary_conditions (double t, double beta[],
double gamma[], double full_u[],
double grid_u[], double dudx[],
int npdes, int ngrids, int left,
int *ires);
#define MIN(X,Y) (X<Y)?X:Y
#define NPDE 2
#define N1 10
#define N2 51
#define N (N1+N2)
#define U(I_,J_) u[I_ * ngrids + J_]
FILE *file1;
int main ()
{
char *state;
int i, j;
int nframes;
double u[(NPDE + 1) * N];
double t0 = 0.0, tout;
double delta_t = 1e-1, tend = 5.0;
int npdes = NPDE, ngrids = N;
double xl = 0.0, xr = 25.0;
double tau = 1.0e-3;
double atol = 1e-2;
double rtol = 0.0;
file1 = fopen ("pde_ex02.out", "w");
imsl_output_file (IMSL_SET_OUTPUT_FILE, file1, 0);
nframes = (int) ((tend + delta_t) / delta_t);
fprintf (file1, " %d\t%d\t%d", npdes, ngrids, nframes);
fprintf (file1, "\t%f\t%f\t%f\t%f\n", xl, xr, t0, tend);
imsl_d_pde_1d_mg_mgr (IMSL_PDE_INITIALIZE, &state,
IMSL_TIME_SMOOTHING, tau,
IMSL_INITIAL_CONDITIONS, initial_conditions, 0);
t0 = 0.0;
tout = delta_t;
do
{
imsl_d_pde_1d_mg (npdes, ngrids, &t0, tout, u, xl,
xr, state, pde_systems, boundary_conditions,
IMSL_RELATIVE_TOLERANCE, rtol,
IMSL_ABSOLUTE_TOLERANCE, atol, 0);
t0 = tout;
fprintf (file1, "%f\n", tout);
for (i = 0; i < npdes + 1; i++)
{
for (j = 0; j < ngrids; j++)
{
fprintf (file1, "%16.10f ", U (i, j));
if (((j + 1) % 4) == 0)
fprintf (file1, "\n");
}
fprintf (file1, "\n");
}
tout = tout + delta_t;
tout = MIN (tout, tend);
}
while (t0 < tend);
imsl_d_pde_1d_mg_mgr (IMSL_PDE_RESET, &state, 0);
fclose (file1);
#undef MIN
#undef NPDE
#undef NFRAMES
#undef N
#undef U
}
static void
initial_conditions (int npdes, int ngrids, double u[])
{
#define U(I_,J_) u[I_* ngrids + J_]
int i, j, i_, n1 = 10, n2 = 51, n;
double dx1, dx2;
double xl = 0.0, xr = 25.0;
n = n1 + n2;
for (i = 0; i < ngrids; i++)
{
U (0, i) = 1.0;
U (1, i) = 0.0;
U (2, i) = 0.0;
}
dx1 = xr / n2;
dx2 = dx1 / n1;
/* grid */
for (i = 1; i <= n1; i++)
{
i_ = i - 1;
U (2, i_) = (i - 1) * dx2;
}
for (i = n1 + 1; i <= n; i++)
{
i_ = i - 1;
U (2, i_) = (i - n1) * dx1;
}
for (i = 0; i < npdes + 1; i++)
{
for (j = 0; j < ngrids; j++)
{
fprintf (file1, "%16.10f ", U (i, j));
if (((j + 1) % 4) == 0)
fprintf (file1, "\n");
}
fprintf (file1, "\n");
}
#undef U
}
static void
pde_systems (double t, double x, int npdes, int ngrids,
double full_u[], double grid_u[], double dudx[],
double *c, double q[], double r[], int *ires)
{
#define C(I_,J_) c[I_ * npdes + J_]
double z;
C (0, 0) = 1.0;
C (1, 0) = 0.0;
C (0, 1) = grid_u[0];
C (1, 1) = 0.0;
r[0] = -grid_u[1];
r[1] = dudx[0];
q[0] = 0.0;
q[1] = grid_u[1] * dudx[0];
return;
#undef C
}
static void
boundary_conditions (double t, double beta[], double gamma[],
double full_u[], double grid_u[], double dudx[],
int npdes, int ngrids, int left, int *ires)
{
double dif;
beta[0] = 0.0;
beta[1] = 0.0;
if (left)
{
dif = exp (-20.0 * t);
gamma[0] = grid_u[0] - dif;
gamma[1] = grid_u[1];
}
else
{
gamma[0] = grid_u[0] - 1.0;
gamma[1] = dudx[1];
}
return;
}
Example 3 - Population Dynamics
This example is from Pennington and Berzins (1994). The system is
This is a notable problem because it involves the unknown
across the entire domain. The software can solve the problem by introducing two dependent algebraic equations:
This leads to the modified system
In the interface to the evaluation of the differential equation and boundary conditions, it is necessary to evaluate the integrals, which are computed with the values of 
on the grid. The integrals are approximated using the trapezoid rule, commensurate with the truncation error in the integrator.
This is a non-linear integro-differential problem involving non-local conditions for the differential equation and boundary conditions. Access to evaluation of these conditions is provided using the optional arguments IMSL_PDE_SYS_W_DATA and IMSL_BOUNDARY_COND_W_DATA. Optional changes are made to use an absolute error tolerance and non-zero time-smoothing. The time-smoothing value 
prevents grid lines from crossing.
#include <stdio.h>
#include <math.h>
#include <imsl.h>
/* prototypes */
static void initial_conditions (int npdes, int ngrids, double u[]);
static void pde_systems (double t, double x, int npdes, int ngrids,
double full_u[], double grid_u[],
double dudx[], double *c, double q[],
double r[], int *ires);
static void boundary_conditions (double t, double beta[],
double gamma[], double full_u[],
double grid_u[], double dudx[],
int npdes, int ngrids, int left,
int *ires);
static double fcn_g (double z, double t);
#define MIN(X,Y) (X<Y)?X:Y
#define NPDE 1
#define N 101
#define U(I_,J_) u[I_ * ngrids + J_]
FILE *file1;
int main ()
{
int i, j, nframes;
double u[(NPDE + 1) * N], mid[N - 1];
int npdes = NPDE, ngrids = N;
double t0 = 0.0, tout;
double delta_t = 1e-1, tend = 5.0, a = 5.0;
char *state;
double xl = 0.0, xr = 5.0;
double *ptr_u;
double tau = 1.0;
double atol = 1e-2;
double rtol = 0.0;
file1 = fopen ("pde_ex03.out", "w");
imsl_output_file (IMSL_SET_OUTPUT_FILE, file1, 0);
nframes = (int) (tend + delta_t) / delta_t;
fprintf (file1, " %d\t%d\t%d", npdes, ngrids, nframes);
fprintf (file1, "\t%f\t%f\t%f\t%f\n", xl, xr, t0, tend);
ptr_u = u;
imsl_d_pde_1d_mg_mgr (IMSL_PDE_INITIALIZE, &state,
IMSL_TIME_SMOOTHING, tau,
IMSL_INITIAL_CONDITIONS, initial_conditions, 0);
tout = delta_t;
fprintf (file1, "%f\n", t0);
do
{
imsl_d_pde_1d_mg (npdes, ngrids, &t0, tout, u, xl,
xr, state, pde_systems, boundary_conditions,
IMSL_RELATIVE_TOLERANCE, rtol,
IMSL_ABSOLUTE_TOLERANCE, atol, 0);
t0 = tout;
if (t0 <= tend)
{
fprintf (file1, "%f\n", tout);
for (i = 0; i < npdes + 1; i++)
{
for (j = 0; j < ngrids; j++)
{
fprintf (file1, "%16.10f ", U (i, j));
if (((j + 1) % 4) == 0)
fprintf (file1, "\n");
}
fprintf (file1, "\n");
}
}
tout = MIN (tout + delta_t, tend);
}
while (t0 < tend);
imsl_d_pde_1d_mg_mgr (IMSL_PDE_RESET, &state, 0);
fclose (file1);
#undef MIN
#undef NPDE
#undef N
#undef XL
#undef XR
#undef U
}
static void
initial_conditions (int npdes, int ngrids, double u[])
{
#define U(I_,J_) u[I_ * ngrids + J_]
#define XL 0.0
#define XR 5.0
int i, j;
double dx, xi;
dx = (XR - XL) / (ngrids - 1);
for (i = 0; i < ngrids; i++)
{
U (0, i) = exp (-U (1, i)) / (2.0 - exp (-XR));
}
for (i = 0; i < npdes + 1; i++)
{
for (j = 0; j < ngrids; j++)
{
fprintf (file1, "%16.10f ", U (i, j));
if (((j + 1) % 4) == 0)
fprintf (file1, "\n");
}
fprintf (file1, "\n");
}
#undef U
#undef XL
#undef XR
}
static void
pde_systems (double t, double x, int npdes, int ngrids,
double full_u[], double grid_u[], double dudx[],
double *c, double q[], double r[], int *ires)
{
#define U(I_,J_) full_u[I_ * ngrids + J_]
double v1;
double sum = 0.0;
int i;
c[0] = 1.0;
r[0] = -1 * grid_u[0];
for (i = 0; i < ngrids - 1; i++)
{
sum += (U (0, i) + U (0, i + 1)) * (U (1, i + 1) - U (1, i));
}
v1 = 0.5 * sum;
q[0] = v1 * grid_u[0];
return;
#undef U
}
static void
boundary_conditions (double t, double beta[], double gamma[],
double full_u[], double grid_u[], double dudx[],
int npdes, int ngrids, int left, int *ires)
{
#define U(I_,J_) full_u[I_ * ngrids + J_]
double v1, v2, mid;
double sum = 0.0;
double sum1 = 0.0, sum2 = 0.0, sum3 = 0.0, sum4 = 0.0;
int i;
for (i = 0; i < ngrids - 1; i++)
{
sum += (U (0, i) + U (0, i + 1)) * (U (1, i + 1) - U (1, i));
mid = 0.5 * (U (1, i) + U (1, i + 1));
sum1 += mid * exp (-mid) *
((U (0, i) + U (0, i + 1)) * (U (1, i + 1) - U (1, i)));
}
if (left)
{
v1 = 0.5 * sum;
v2 = 0.5 * sum1;
beta[0] = 0.0;
gamma[0] = fcn_g (1.0, t) * v1 * v2 /
((v1 + 1.0) * (v1 + 1.0)) - grid_u[0];
}
else
{
beta[0] = 0.0;
gamma[0] = dudx[0];
}
return;
#undef U
}
static double
fcn_g (double z, double t)
{
double g, a = 5.0;
g = 4.0 * z * (2.0 - 2.0 * exp (-a) + exp (-t)) *
(2.0 - 2.0 * exp (-a) + exp (-t));
g = g / ((1.0 - exp (-a)) * (1.0 - (1.0 + 2.0 * a) *
exp (-2.0 * a)) * (1.0 - exp (-a) + exp (-t)));
return g;
}
Example 4 - A Model in Cylindrical Coordinates
This example is from Blom and Zegeling (1994). The system models a reactor-diffusion problem:
The axial direction 
is treated as a time coordinate. The radius 
is treated as the single space variable.
This is a non-linear problem in cylindrical coordinates. Our example illustrates assigning 
in Equation 2. We provide the optional argument IMSL_CYL_COORDINATES that resets this value from its default, 
.
#include <stdio.h>
#include <math.h>
#include <imsl.h>
/* prototypes */
static void initial_conditions (int npdes, int ngrids, double t[]);
static void pde_systems (double t, double x, int npdes, int ngrids,
double u[], double grid_u[], double dudx[], double *c,
double q[], double r[], int *ires);
static void boundary_conditions (double t, double beta[],
double gamma[], double u[], double grid_u[], double dudx[],
int npdes, int ngrids, int left, int *ires);
#define MIN(X,Y) (X<Y)?X:Y
#define NPDE 1
#define N 41
#define T(I_,J_) t[I_ * ngrids + J_]
int main ()
{
int i, j, ido;
int nframes;
double t[(NPDE + 1) * N];
double z0 = 0.0, zout;
double dx1, dx2, diff;
double delta_z = 1e-1, zend = 1.0, zmax = 1.0;
double beta = 1e-4, gamma = 1.0, eps = 1e-1;
char *state;
int npdes = NPDE, ngrids = N;
double xl = 0.0, xr = 1.0;
FILE *file1;
int m = 1;
file1 = fopen ("pde_ex04.out", "w");
imsl_output_file (IMSL_SET_OUTPUT_FILE, file1, 0);
nframes = (int) ((zend + delta_z) / delta_z) - 1;
fprintf (file1, " %d\t%d\t%d", npdes, N, nframes);
fprintf (file1, "\t%f\t%f\t%f\t%f\n", xl, xr, z0, zend);
imsl_d_pde_1d_mg_mgr (IMSL_PDE_INITIALIZE, &state, IMSL_CYL_COORDINATES, 0);
initial_conditions (npdes, ngrids, t);
zout = delta_z;
do
{
imsl_d_pde_1d_mg (npdes, ngrids, &z0, zout, t, xl,
xr, state, pde_systems, boundary_conditions, 0);
z0 = zout;
if (z0 <= zend)
{
fprintf (file1, "%f\n", zout);
for (i = 0; i < npdes + 1; i++)
{
for (j = 0; j < ngrids; j++)
{
fprintf (file1, "%16.10f ", T (i, j));
if (((j + 1) % 4) == 0)
fprintf (file1, "\n");
}
fprintf (file1, "\n");
}
}
zout = MIN ((zout + delta_z), zend);
}
while (z0 < zend);
imsl_d_pde_1d_mg_mgr (IMSL_PDE_RESET, &state, 0);
fclose (file1);
#undef MIN
#undef NPDE
#undef N
#undef T
}
static void
initial_conditions (int npdes, int ngrids, double t[])
{
#define T(I_,J_) t[I_ * ngrids + J_]
int i;
for (i = 0; i < ngrids; i++)
{
T (0, i) = 0.0;
}
#undef T
}
static void
pde_systems (double t, double x, int npdes, int ngrids, double u[],
double grid_u[], double dudx[], double *c,
double q[], double r[], int *ires)
{
#define C(I_,J_) c[I_ * npdes + J_]
static double beta = 01e-4, gamma = 1.0, eps = 1e-1;
C (0, 0) = 1.0;
r[0] = beta * dudx[0];
q[0] = -1.0 * gamma * exp (grid_u[0] / (1.0 + eps * grid_u[0]));
return;
#undef C
}
static void
boundary_conditions (double t, double beta[], double gamma[],
double u[], double grid_u[], double dudx[],
int npdes, int ngrids, int left, int *ires)
{
if (left)
{
beta[0] = 1.0;
gamma[0] = 0.0;
}
else
{
beta[0] = 0.0;
gamma[0] = grid_u[0];
}
return;
}
Example 5 - A Flame Propagation Model
This example is presented more fully in Verwer, et al., (1989). The system is a normalized problem relating mass density 
and temperature 
:
This is a non-linear problem. The example shows the model steps for replacing the banded solver in the software with one of the user’s choice. Following the computation of the matrix factorization in imsl_lin_sol_gen_band (see Chapter 1, Linear Systems), we declare the system to be singular when the reciprocal of the condition number is smaller than the working precision. This choice is not suitable for all problems. Attention must be given to detecting a singularity when this option is used.
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
#include <imsl.h>
/* prototypes */
static void initial_conditions (int npdes, int ngrids, double u[]);
static void pde_systems (double t, double x, int npdes, int ngrids,
double u[], double grid_u[], double dudx[],
double *c, double q[], double r[], int *ires);
static void boundary_conditions (double t, double beta[],
double gamma[], double u[],
double grid_u[], double dudx[],
int npdes, int ngrids, int left,
int *ires);
static int fac (int neq, int iband, double *a);
static void sol (int neq, int iband, double *g, double *y);
static double fcn (double z);
int *ipvt = NULL;
double *factor = NULL;
#define MIN(X,Y) (X<Y)?X:Y
#define NPDE 2
#define N 40
#define NEQ ((NPDE+1)*N)
#define U(I_,J_) u[I_ * ngrids + J_]
int main ()
{
int i, j, nframes;
double u[(NPDE + 1) * N];
double t0 = 0.0, tout;
double delta_t = 1e-4, tend = 6e-3;
char *state;
int npdes = NPDE, ngrids = N;
double xl = 0.0, xr = 1.0;
FILE *file1;
double work[NEQ], rcond;
double xmax = 1.0, beta = 4.0, gamma = 3.52e6;
int max_bdf_order = 5;
file1 = fopen ("pde_ex05.out", "w");
imsl_output_file (IMSL_SET_OUTPUT_FILE, file1, 0);
nframes = (int) ((tend + delta_t) / delta_t) - 1;
fprintf (file1, " %d\t%d\t%d", npdes, ngrids, nframes);
fprintf (file1, "\t%f\t%f\t%f\t%f\n", xl, xr, t0, tend);
initial_conditions (npdes, ngrids, u);
imsl_d_pde_1d_mg_mgr (IMSL_PDE_INITIALIZE, &state,
IMSL_MAX_BDF_ORDER, max_bdf_order,
IMSL_USER_FACTOR_SOLVE, fac, sol, 0);
tout = delta_t;
do
{
imsl_d_pde_1d_mg (npdes, ngrids, &t0, tout, u, xl,
xr, state, pde_systems, boundary_conditions, 0);
t0 = tout;
if (t0 <= tend)
{
fprintf (file1, "%f\n", tout);
for (i = 0; i < npdes + 1; i++)
{
for (j = 0; j < ngrids; j++)
{
fprintf (file1, "%16.10f ", U (i, j));
if (((j + 1) % 4) == 0)
fprintf (file1, "\n");
}
}
}
tout = MIN ((tout + delta_t), tend);
}
while (t0 < tend);
imsl_d_pde_1d_mg_mgr (IMSL_PDE_RESET, &state, 0);
fclose (file1);
if (factor != NULL)
{
imsl_free (factor);
}
if (ipvt != NULL)
{
imsl_free (ipvt);
}
}
#undef MIN
#undef U
static void
initial_conditions (int npdes, int ngrids, double u[])
{
#define U(I_,J_) u[I_ * ngrids + J_]
int i;
for (i = 0; i < ngrids; i++)
{
U (0, i) = 1.0;
U (1, i) = 2e-1;
}
#undef U
}
static void
pde_systems (double t, double x, int npdes, int ngrids, double u[],
double grid_u[], double dudx[], double *c,
double q[], double r[], int *ires)
{
#define C(I_,J_) c[I_ * npdes + J_]
C (0, 0) = 1.0;
C (0, 1) = 0.0;
C (1, 0) = 0.0;
C (1, 1) = 1.0;
r[0] = dudx[0];
r[1] = dudx[1];
q[0] = grid_u[0] * fcn (grid_u[1]);
q[1] = -1.0 * q[0];
return;
#undef C
}
static void
boundary_conditions (double t, double beta[], double gamma[],
double u[], double grid_u[], double dudx[],
int npdes, int ngrids, int left, int *ires)
{
if (left)
{
beta[0] = 0.0;
beta[1] = 0.0;
gamma[0] = dudx[0];
gamma[1] = dudx[1];
}
else
{
beta[0] = 1.0;
gamma[0] = 0.0;
beta[1] = 0.0;
if (t >= 2e-4)
{
gamma[1] = 12e-1;
}
else
{
gamma[1] = 2e-1 + 5e3 * t;
}
gamma[1] -= grid_u[1];
}
return;
}
/* Factor the banded matrix. This is the same solver used
* internally but that is not required. A user can substitute
* one of their own.
* Note: Allowing lin_sol_gen_band to allocate ipvt and factor
* variables, then use in sol function.
*/
static int
fac (int neq, int iband, double *a)
{
double rcond, panic_flag;
int i, j;
double b[NEQ];
/* Free factor and pivot sequence if previously allocated. */
if (factor != NULL)
{
imsl_free (factor);
factor = NULL;
}
if (ipvt != NULL)
{
imsl_free (ipvt);
ipvt = NULL;
}
imsl_d_lin_sol_gen_band (neq, a, iband, iband, b,
IMSL_FACTOR, &ipvt, &factor,
IMSL_FACTOR_ONLY, IMSL_CONDITION, &rcond, 0);
panic_flag = 0;
if (1.0 / rcond <= imsl_d_machine (4))
panic_flag = 3;
return panic_flag;
}
static void
sol (int neq, int iband, double *g, double *y)
{
imsl_d_lin_sol_gen_band (neq, (double *) NULL, iband, iband, g,
IMSL_SOLVE_ONLY,
IMSL_FACTOR_USER, ipvt, factor,
IMSL_RETURN_USER, y, 0);
return;
}
static double
fcn (double z)
{
double beta = 4.0, gamma = 3.52e6;
return gamma * exp (-1.0 * beta / z);
}
Example 6 - A ‘Hot Spot’ Model
This example is presented more fully in Verwer, et al., (1989). The system is a normalized problem relating the temperature 
, of a reactant in a chemical system. The formula for 
is equivalent to their example.
This is a non-linear problem. The output shows a case where a rapidly changing front, or hot-spot, develops after a considerable way into the integration. This causes rapid change to the grid. An option sets the maximum order BDF formula from its default value of 2 to the theoretical stable maximum value of 5.
#include <stdio.h>
#include <math.h>
#include <imsl.h>
/* prototypes */
static void initial_conditions (int npdes, int ngrids, double u[]);
static void pde_systems (double t, double x, int npdes, int ngrids,
double full_u[], double grid_u[],
double dudx[], double *c, double q[],
double r[], int *ires);
static void boundary_conditions (double t, double beta[],
double gamma[], double full_u[],
double grid_u[], double dudx[],
int npdes, int ngrids, int left,
int *ires);
static double fcn_h (double z);
#define MIN(X,Y) (X<Y)?X:Y
#define NPDE 1
#define N 80
#define U(I_,J_) u[I_ * ngrids + J_]
int main ()
{
int i, j, nframes;
double u[(NPDE + 1) * N];
double t0 = 0.0, tout;
double delta_t = 1e-2, tend = 29e-2;
double u0 = 1.0, u1 = 0.0, tdelta = 1e-1, tol = 29e-2;
double a = 1.0, delta = 20.0, r = 5.0;
char *state;
int npdes = NPDE, ngrids = N;
double xl = 0.0, xr = 1.0;
FILE *file1;
int max_bdf_order = 5;
file1 = fopen ("pde_ex06.out", "w");
imsl_output_file (IMSL_SET_OUTPUT_FILE, file1, 0);
nframes = (int) ((tend + delta_t) / delta_t) - 1;
fprintf (file1, " %d\t%d\t%d", npdes, ngrids, nframes);
fprintf (file1, "\t%f\t%f\t%f\t%f\n", xl, xr, t0, tend);
initial_conditions (npdes, ngrids, u);
imsl_d_pde_1d_mg_mgr (IMSL_PDE_INITIALIZE, &state,
IMSL_MAX_BDF_ORDER, max_bdf_order, 0);
tout = delta_t;
do
{
imsl_d_pde_1d_mg (npdes, ngrids, &t0, tout, u, xl,
xr, state, pde_systems, boundary_conditions, 0);
t0 = tout;
if (t0 <= tend)
{
fprintf (file1, "%f\n", tout);
for (i = 0; i < npdes + 1; i++)
{
for (j = 0; j < ngrids; j++)
{
fprintf (file1, "%16.10f ", U (i, j));
if (((j + 1) % 4) == 0) fprintf (file1, "\n");
}
}
}
tout = MIN ((tout + delta_t), tend);
}
while (t0 < tend);
imsl_d_pde_1d_mg_mgr (IMSL_PDE_RESET, &state, 0);
fclose (file1);
#undef MIN
#undef NPDE
#undef N
#undef U
}
static void
initial_conditions (int npdes, int ngrids, double u[])
{
#define U(I_,J_) u[I_ * ngrids + J_]
int i;
for (i = 0; i < ngrids; i++)
{
U (0, i) = 1.0;
}
#undef U
}
static void
pde_systems (double t, double x, int npdes, int ngrids,
double full_u[], double grid_u[], double dudx[],
double *c, double q[], double r[], int *ires)
{
#define C(I_,J_) c[I_ * npdes + J_]
c[0] = 1.0;
r[0] = dudx[0];
q[0] = -fcn_h (grid_u[0]);
return;
#undef C
}
static void
boundary_conditions (double t, double beta[], double gamma[],
double full_u[], double grid_u[], double dudx[],
int npdes, int ngrids, int left, int *ires)
{
if (left)
{
beta[0] = 0.0;
gamma[0] = dudx[0];
}
else
{
beta[0] = 0.0;
gamma[0] = grid_u[0] - 1.0;
}
return;
}
static double
fcn_h (double z)
{
double a = 1.0, delta = 2e1, r = 5.0;
return (r / (a * delta)) * (1.0 + a - z) *
exp (-delta * (1.0 / z - 1.0));
}
This example is presented more fully in Verwer, et al., (1989). The system is a normalized problem relating the interaction of two waves, 
and 
moving in opposite directions. The waves meet and reduce in amplitude, due to the non-linear terms in the equation. Then they separate and travel onward, with reduced amplitude.
This is a non-linear system of first order equations.
#include <stdio.h>
#include <math.h>
#include <imsl.h>
/* prototypes */
static void initial_conditions (int npdes, int ngrids, double u[]);
static void pde_systems (double t, double x, int npdes, int ngrids,
double full_u[], double grid_u[],
double dudx[], double *c, double q[],
double r[], int *ires);
static void boundary_conditions (double t, double beta[],
double gamma[], double full_u[],
double grid_u[], double dudx[],
int npdes, int ngrids, int left,
int *ires);
#define MIN(X,Y) (X<Y)?X:Y
#define NPDE 2
#define N 50
#define XL (-0.5)
#define XR 0.5
#define U(I_,J_) u[I_ * ngrids + J_]
FILE *file1;
int main ()
{
int i, j, nframes;
double u[(NPDE + 1) * N];
double t0 = 0.0, tout;
double delta_t = 5e-2, tend = 5e-1;
char *state;
int npdes = NPDE, ngrids = N;
double xl = XL, xr = XR;
double tau = 1e-3;
double atol = 1e-3;
double rtol = 0.0;
int max_bdf_order = 3;
file1 = fopen ("pde_ex07.out", "w");
imsl_output_file (IMSL_SET_OUTPUT_FILE, file1, 0);
nframes = (int) ((tend + delta_t) / delta_t);
fprintf (file1, " %d\t%d\t%d", npdes, ngrids, nframes);
fprintf (file1, "\t%f\t%f\t%f\t%f\n", xl, xr, t0, tend);
imsl_d_pde_1d_mg_mgr (IMSL_PDE_INITIALIZE, &state,
IMSL_TIME_SMOOTHING, tau,
IMSL_MAX_BDF_ORDER, max_bdf_order,
IMSL_INITIAL_CONDITIONS, initial_conditions, 0);
fprintf (file1, "%f\n", t0);
tout = delta_t;
do
{
imsl_d_pde_1d_mg (npdes, ngrids, &t0, tout, u, xl,
xr, state, pde_systems, boundary_conditions,
IMSL_RELATIVE_TOLERANCE, rtol,
IMSL_ABSOLUTE_TOLERANCE, atol, 0);
t0 = tout;
if (t0 <= tend)
{
fprintf (file1, "%f\n", tout);
for (i = 0; i < npdes + 1; i++)
{
for (j = 0; j < ngrids; j++)
{
fprintf (file1, "%16.10f ", U (i, j));
if (((j + 1) % 4) == 0)
fprintf (file1, "\n");
}
fprintf (file1, "\n");
}
}
tout = MIN ((tout + delta_t), tend);
}
while (t0 < tend);
imsl_d_pde_1d_mg_mgr (IMSL_PDE_RESET, &state, 0);
fclose (file1);
#undef MIN
#undef NPDE
#undef N
#undef XL
#undef XR
#undef U
}
static void
initial_conditions (int npdes, int ngrids, double u[])
{
#define U(I_,J_) u[I_ * ngrids + J_]
#define XL -0.5
#define XR 0.5
int i, j;
double _pi, pulse;
double dx, xi;
_pi = imsl_d_constant("pi",0);
for (i = 0; i < ngrids; i++)
{
pulse = (0.5 * (1.0 + cos (10.0 * _pi * U (npdes, i))));
U (0, i) = pulse;
U (1, i) = pulse;
}
for (i = 0; i < ngrids; i++)
{
if ((U (npdes, i) < -3e-1) || (U (npdes, i) > -1e-1))
{
U (0, i) = 0.0;
}
if ((U (npdes, i) < 1e-1 || U (npdes, i) > 3e-1))
{
U (1, i) = 0.0;
}
}
for (i = 0; i < npdes + 1; i++)
{
for (j = 0; j < ngrids; j++)
{
fprintf (file1, "%16.10f ", U (i, j));
if (((j + 1) % 4) == 0)
fprintf (file1, "\n");
}
fprintf (file1, "\n");
}
#undef XL
#undef XR
#undef U
}
static void
pde_systems (double t, double x, int npdes, int ngrids,
double full_u[], double grid_u[], double dudx[],
double *c, double q[], double r[], int *ires)
{
#define C(I_,J_) c[I_ * npdes + J_]
C (0, 0) = 1.0;
C (0, 1) = 0.0;
C (1, 0) = 0.0;
C (1, 1) = 1.0;
r[0] = -1.0 * grid_u[0];
r[1] = grid_u[1];
q[0] = 100.0 * grid_u[0] * grid_u[1];
q[1] = q[0];
return;
#undef C
}
static void
boundary_conditions (double t, double beta[], double gamma[],
double full_u[], double grid_u[], double dudx[],
int npdes, int ngrids, int left, int *ires)
{
beta[0] = 0.0;
beta[1] = 0.0;
gamma[0] = grid_u[0];
gamma[1] = grid_u[1];
return;
}
The value of a European “call option,” 
, with exercise price 
and expiration date 
, satisfies the “asset-or-nothing payoff ” 
. Prior to expiration 
is estimated by the Black-Scholes differential equation 
. The parameters in the model are the risk-free interest rate, 
, and the stock volatility, 
. The boundary conditions are 
and 
. This development is described in Wilmott, et al. (1995), pages 41-57. There are explicit solutions for this equation based on the Normal Curve of Probability. The normal curve, and the solution itself, can be efficiently computed with the IMSL function imsl_f_normal_cdf, see Chapter 9, “Special Functions.” With numerical integration the equation itself or the payoff can be readily changed to include other formulas, 
, and corresponding boundary conditions. We use 
.
This is a linear problem but with initial conditions that are discontinuous. It is necessary to use a positive time-smoothing value to prevent grid lines from crossing. We have used an absolute tolerance of 
. In $US, this is one-tenth of a cent.
#include <stdio.h>
#include <imsl.h>
/* prototypes */
static void initial_conditions (int npdes, int ngrids, double u[]);
static void pde_systems (double t, double x, int npdes, int ngrids,
double full_u[], double grid_u[], double dudx[], double *c,
double q[], double r[], int *ires);
static void boundary_conditions (double t, double beta[], double gamma[],
double full_u[], double grid_u[], double dudx[], int npdes,
int ngrids, int left, int *ires);
#define MIN(X,Y) (X<Y)?X:Y
#define NPDE 1
#define N 100
#define XL 0.0
#define XR 150.0
#define U(I_,J_) u[I_ * ngrids + J_]
int main ()
{
int i, j, nframes;
double u[(NPDE + 1) * N];
double t0 = 0.0, tout, xval;
double delta_t = 25e-3, tend = 25e-2;
double xmax = 150.0;
char *state;
int npdes = NPDE, ngrids = N;
double xl = XL, xr = XR;
FILE *file1;
double tau = 5e-3;
double atol = 1e-2;
double rtol = 0.0;
int max_bdf_order = 5;
file1 = fopen ("pde_ex08.out", "w");
imsl_output_file (IMSL_SET_OUTPUT_FILE, file1, 0);
nframes = (int) ((tend + delta_t) / delta_t);
fprintf (file1, " %d\t%d\t%d", npdes, ngrids, nframes);
fprintf (file1, "\t%f\t%f\t%f\t%f\n", xl, xr, t0, tend);
initial_conditions (npdes, ngrids, u);
imsl_d_pde_1d_mg_mgr (IMSL_PDE_INITIALIZE, &state,
IMSL_TIME_SMOOTHING, tau,
IMSL_MAX_BDF_ORDER, max_bdf_order, 0);
tout = delta_t;
do
{
imsl_d_pde_1d_mg (npdes, ngrids, &t0, tout, u, xl,
xr, state, pde_systems, boundary_conditions,
IMSL_RELATIVE_TOLERANCE, rtol,
IMSL_ABSOLUTE_TOLERANCE, atol, 0);
t0 = tout;
if (t0 <= tend)
{
fprintf (file1, "%f\n", tout);
for (i = 0; i < npdes + 1; i++)
{
for (j = 0; j < ngrids; j++)
{
fprintf (file1, "%16.10f ", U (i, j));
if (((j + 1) % 4) == 0)
fprintf (file1, "\n");
}
}
}
tout = MIN ((tout + delta_t), tend);
}
while (t0 < tend);
imsl_d_pde_1d_mg_mgr (IMSL_PDE_RESET, &state, 0);
fclose (file1);
#undef MIN
#undef NPDE
#undef N
#undef XL
#undef XR
#undef U
}
static void
initial_conditions (int npdes, int ngrids, double u[])
{
#define U(I_,J_) u[I_ * ngrids + J_]
#define XL 0.0
#define XR 150.0
int i;
double dx, xi, xval, e = 100.0;
dx = (XR - XL) / (ngrids - 1);
for (i = 0; i < ngrids; i++)
{
xi = XL + i * dx;
if (xi <= e)
{
U (0, i) = 0.0;
}
else
{
U (0, i) = xi;
}
}
#undef U
#undef XL
#undef XR
}
static void
pde_systems (double t, double x, int npdes, int ngrids,
double full_u[], double grid_u[], double dudx[], double *c,
double q[], double r[], int *ires)
{
double sigsq, sigma = 2e-1, rr = 8e-2;
sigsq = sigma * sigma;
c[0] = 1.0;
r[0] = dudx[0] * x * x * sigsq * 0.5;
q[0] = -(rr - sigsq) * x * dudx[0] + rr * grid_u[0];
return;
}
static void
boundary_conditions (double t, double beta[], double gamma[],
double full_u[], double grid_u[], double dudx[],
int npdes, int ngrids, int left, int *ires)
{
if (left)
{
beta[0] = 0.0;
gamma[0] = grid_u[0];
}
else
{
beta[0] = 0.0;
gamma[0] = dudx[0] - 1.0;
}
return;
}
PRO PDE_1d_mg_plot, FILENAME = filename, PAUSE = pause
;
if keyword_set(FILENAME) then file = filename else file = "res.dat"
if keyword_set(PAUSE) then twait = pause else twait = .1
;
; Define floating point variables that will be read
; from the first line of the data file.
xl = 0D0
xr = 0D0
t0 = 0D0
tlast = 0D0
;
; Open the data file and read in the problem parameters.
openr, lun, filename, /get_lun
readf, lun, npde, np, nt, xl, xr, t0, tlast
; Define the arrays for the solutions and grid.
u = dblarr(nt, npde, np)
g = dblarr(nt, np)
times = dblarr(nt)
;
; Define a temporary array for reading in the data.
tmp = dblarr(np)
t_tmp = 0D0
;
; Read in the data.
for i = 0, nt-1 do begin ; For each step in time
readf, lun, t_tmp
times(i) = t_tmp
for k = 0, npde-1 do begin ; For each PDE:
rmf, lun, tmp
u(i,k,*) = tmp ; Read in the components.
end
rmf, lun, tmp
g(i,*) = tmp ; Read in the grid.
end
;
; Close the data file and free the unit.
close, lun
free_lun, lun
;
; We now have all of the solutions and grids.
;
; Delete any window that is currently open.
while (!d.window NE -1) do WDELETE
;
; Open two windows for plotting the solutions
; and grid.
window, 0, xsize = 550, ysize = 420
window, 1, xsize = 550, ysize = 420
;
; Plot the grid.
wset, 0
plot, [xl, xr], [t0, tlast], /nodata, ystyle = 1, $
title = "Grid Points", xtitle = "X", ytitle = "Time"
for i = 0, np-1 do begin
oplot, g(*, i), times, psym = -1
end
;
; Plot the solution(s):
wset, 1
for k = 0, npde-1 do begin
umin = min(u(*,k,*))
umax = max(u(*,k,*))
for i = 0, nt-1 do begin
title = strcompress("U_"+string(k+1), /remove_all)+ $
" at time "+string(times(i))
plot, g(i, *), u(i,k,*), ystyle = 1, $
title = title, xtitle = "X", $
ytitle = strcompress("U_"+string(k+1), /remove_all), $
xr = [xl, xr], yr = [umin, umax], $
psym = -4
wait, twait
end
end
end
Fatal Errors
|
IMSL_STOP_USER_FCN |
Request from user supplied function to stop algorithm. |
