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
inttask (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(floatt, floatx, intnpdes, intngrids, 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 (floatt, float*beta, float*gamma, float *full_u, float *grid_u, float*dudx, intnpdes, intgrids, intleft, 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
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, floattau, (Input) Resets the value of the parameter , described above. Default: .
IMSL_SPATIAL_SMOOTHING, floatkappa, (Input) Resets the value of the parameter , described above. Default: .
IMSL_MONITOR_REGULARIZING, floatalpha, (Input) Resets the value of the parameter , described above. Default: .
IMSL_MAX_BDF_ORDER, intmax_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(intneq, 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(intneq,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(intnpdes, intngrids, 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(intnpdes, intngrids, 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
IMSL_RELATIVE_TOLERANCE, floatrtol, (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, floatatol, (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(floatt, floatx, intnpdes, intngrids, 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 (floatt,float*beta, float*gamma, float*full_u, float*grid_u, float*dudx, intnpdes, intngrids, intleft,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
Remarks on the 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.
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,
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,
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,
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,
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,
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,
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,
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,
