Solves a system of partial differential equations of the form ut = f(x, t, u, ux, uxx) using the method of lines. The solution is represented with cubic Hermite polynomials.
Note: imsl_f_modified_method_of_lines replaces deprecated function imsl_f_pde_method_of_lines.
Synopsis
#include <imsl.h>
void imsl_f_modified_method_of_lines_mgr (int task, void **state, ..., 0)
void imsl_f_modified_method_of_lines (int npdes, float *t, float tend, int nx, float xbreak[], float y[], void *state, void fcn_ut(), void fcn_bc())
The type double functions are imsl_d_modified_method_of_lines_mgr and imsl_d_modified_method_of_lines.
Required Arguments for imsl_f_modified_method_of_lines_mgr
int task
(Input)
This function must be called with task set to IMSL_PDE_INITIALIZE to
set up memory and default values prior to solving a problem and with task equal
to IMSL_PDE_RESET to
clean up after it has solved.
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_modified_method_of_lines
int npdes
(Input)
Number of differential equations.
float *t
(Input/Output)
Independent variable. On input, t supplies the initial
time, t0. On output, t is set to the value
to which the integration has been updated. Normally, this new value is tend.
float tend
(Input)
Value of t = tend at which the
solution is desired.
int nx
(Input)
Number of mesh points or lines.
float xbreak[]
(Input)
Array of length nx containing the
breakpoints for the cubic Hermite splines used in the x discretization.
The points in xbreak must be
strictly increasing. The values xbreak[0] and xbreak[nx − 1] are
the endpoints of the interval.
float y[]
(Input/Output)
Array of length npdes by nx containing the
solution. The array y contains the
solution as y[k,i] =
uk(x, tend) at
x = xbreak[i]. On input, y
contains the initial values. It must satisfy the boundary conditions. On output,
y contains the
computed solution.
void *state
(Input/Output)
The current state of the PDE solution is held in a structure
pointed to by state. It must be initialized by a call to imsl_f_modified_method_of_lines_mgr.
It cannot be directly manipulated.
void fcn_ut(int npdes, float x, float t, float u[], float ux[], float uxx[], float ut[])
User-supplied
function to evaluate ut.
int npdes
(Input)
Number of equations.
float x
(Input)
Space variable, x.
float t
(Input)
Time variable, t.
float u[]
(Input)
Array of length npdes containing the
dependent values, u.
float ux[]
(Input)
Array of length npdes containing the
first derivatives, ux.
float uxx[]
(Input)
Array of length npdes containing the
second derivative, uxx.
float ut[]
(Output)
Array of length npdes containing the
computed derivatives ut.
void fcn_bc(int npdes, float x, float t, float alpha[], float beta[],
float gamma[])
User-supplied
function to evaluate the boundary conditions. The boundary conditions accepted
by imsl_f_modified_method_of_lines
are
NOTE: Users must supply the values αk, βk, and γk.
int npdes
(Input)
Number of equations.
float x
(Input)
Space variable, x.
float t
(Input)
Time variable, t.
float alpha[]
(Output)
Array of length npdes containing the
αk
values.
float beta[]
(Output)
Array of length npdes containing the
βk
values.
float gamma[]
(Output)
Array of length npdes containing the
values of γk.
Synopsis with Optional Arguments
#include <imsl.h>
void imsl_f_modified_method_of_lines_mgr (int task, void **state,
IMSL_TOL, float tol,
IMSL_HINIT, float hinit,
IMSL_INITIAL_VALUE_DERIVATIVE, float initial_deriv[],
IMSL_HTRIAL, float *htrial,
IMSL_FCN_UT_W_DATA, void fcn_ut(), void *data,
IMSL_FCN_BC_W_DATA, void fcn_bc(), void *data,
0)
Optional Arguments
IMSL_TOL, float tol
(Input)
Differential equation error tolerance. An attempt is made to control
the local error in such a way that the global relative error is proportional to
tol.
Default:
tol = 100.0*imsl_f_machine(4)
IMSL_HINIT, float hinit
(Input)
Initial step size in the t integration. This value must be
nonnegative. If hinit is zero, an
initial step size of 0.001|tend - t0| will be
arbitrarily used. The step will be applied in the direction of
integration.
Default: hinit = 0.0
IMSL_INITIAL_VALUE_DERIVATIVE, float initial_deriv[]
(Input/Output)
Supply the derivative values ux(x,
t0)
in initial_deriv, an
array of length npdes by nx. This derivative
information is input as
The array initial_deriv contains the derivative values as output:
Default: Derivatives are computed using cubic spline interpolation.
IMSL_HTRIAL, float *htrial
(Output)
Return the current trial step size.
IMSL_FCN_UT_W_DATA, void fcn_ut(int npdes, float x, float t, float u[], float ux[], float uxx[], float ut[], void *data),
void *data
(Input)
User-supplied function to evaluate ut, 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_FCN_BC_W_DATA, void fcn_bc(int npdes, float x, float t, float alpha[], float beta[],
float gamma[], void *data),
void *data
(Input)
User-supplied function to evaluate 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 M = npdes, N = nx and xi = xbreak(I). The function imsl_f_modified_method_of_lines uses the method of lines to solve the partial differential equation system
with the initial conditions
and the boundary conditions
for k = 1, …, M.
Cubic Hermite polynomials are used in the x variable approximation so that the trial solution is expanded in the series
where φi(x) and ψi(x) are the standard basis functions for the cubic Hermite polynomials with the knots x1 < x2 < … < xN. These are piecewise cubic polynomials with continuous first derivatives. At the breakpoints, they satisfy
According to the collocation method, the coefficients of the approximation are obtained so that the trial solution satisfies the differential equation at the two Gaussian points in each subinterval,
for j = 1, …, N. The collocation approximation to the differential equation is
for k = 1, …, M and j = 1, …, 2(N − 1).
This is a system of 2M(N − 1) ordinary differential equations in 2M N unknown coefficient functions, ak,i and bk,i. This system can be written in the matrix-vector form as A dc/dt = F (t, c) with c(t0) = c0 where c is a vector of coefficients of length 2M N and c0 holds the initial values of the coefficients. The last 2M equations are obtained from the boundary conditions.
If αk = βk = 0, it is assumed that no boundary condition is desired for the k-th unknown at the left endpoint. A similar comment holds for the right endpoint. Thus, collocation is done at the endpoint. This is generally a useful feature for systems of first-order partial differential equations.
The input/output array Y contains the values of the ak,i. The initial values of the bk,i are obtained by using the IMSL cubic spline function imsl_f_cub_spline_interp_e_cnd (Chapter 3, “Interpolation and Approximation”) to construct functions
such that
The IMSL function imsl_f_cub_spline_value, Chapter 3, “Interpolation and Approximation;” is used to approximate the values
Optional argument IMSL_INITIAL_VALUE_DERIVATIVE allows the user to provide the initial values of bk,i.
The order of matrix A is 2M N and its maximum bandwidth is 6M - 1. The band structure of the Jacobian of F with respect to c is the same as the band structure of A. This system is solved using a modified version of imsl_f_ode_adams_gear. Numerical Jacobians are used exclusively. Gear’s BDF method is used as the default because the system is typically stiff. For more details, see Sewell (1982).
Four examples of PDEs are now presented that illustrate how users can interface their problems with IMSL PDE solving software. The examples are small and not indicative of the complexities that most practitioners will face in their applications. A set of seven sample application problems, some of them with more than one equation, is given in Sincovec and Madsen (1975). Two further examples are given in Madsen and Sincovec (1979).
Examples
Example 1
The normalized linear diffusion PDE, ut = uxx, 0 ≤ x ≤ 1, t > t0, is solved. The initial values are t0 = 0, u(x, t0) = u0 = 1. There is a “zero-flux” boundary condition at x = 1, namely ux(1, t) = 0, (t > t0). The boundary value of u(0, t) is abruptly changed from u0 to the value 0, for t >0.
When the boundary conditions are discontinuous, or incompatible with the initial conditions such as in this example, it may be important to use double precision.
#include <imsl.h>
#include <stdio.h>
void fcnut(int, float, float, float *, float *, float *, float *);
void fcnbc(int, float, float, float *, float *, float *);
int main() {
int npdes = 1, nx = 8, i, j, nstep = 10;
float hinit, tol, t = 0.0, tend, xbreak[8], y[8];
char title[50];
void *state;
/* Set breakpoints and initial conditions */
for (i = 0; i < nx; i++) {
xbreak[i] = (float) i / (float) (nx - 1);
y[i] = 1.0;
}
/* Initialize the solver */
tol = 10.e-4;
hinit = 0.01 * tol;
imsl_f_modified_method_of_lines_mgr(IMSL_PDE_INITIALIZE, &state,
IMSL_TOL, tol,
IMSL_HINIT, hinit,
0);
for (j = 1; j <= nstep; j++) {
tend = (float) j / (float) nstep;
tend *= tend;
/* Solve the problem */
imsl_f_modified_method_of_lines(npdes, &t, tend, nx, xbreak, y,
state, fcnut, fcnbc);
/* Print results at current t=tend */
sprintf(title, "solution at t = %4.2f", t);
imsl_f_write_matrix(title, npdes, nx, y, 0);
}
}
void fcnut(int npdes, float x, float t, float *u, float *ux,
float *uxx, float *ut) {
/* Define the PDE */
*ut = *uxx;
}
void fcnbc(int npdes, float x, float t, float *alpha, float *beta,
float *gam) {
float delta = 0.09, u0 = 1.0, u1 = 0.1;
/* Define boundary conditions */
if (x == 0.0) {
/* These are for x = 0 */
*alpha = 1.0;
*beta = 0.0;
*gam = u1;
/* If in the boundary layer, compute
nonzero gamma */
if (t <= delta)
*gam = u0 +(u1 - u0)* t/delta;
} else {
/* These are for x = 1 */
*alpha = 0.0;
*beta = 1.0;
*gam = 0.0;
}
}
Output
solution at t = 0.01
1 2 3 4 5 6
0.900 0.985 0.999 1.000 1.000 1.000
7 8
1.000 1.000
solution at t = 0.04
1 2 3 4 5 6
0.600 0.834 0.941 0.982 0.996 0.999
7 8
1.000 1.000
solution at t = 0.09
1 2 3 4 5 6
0.1003 0.4906 0.7304 0.8673 0.9395 0.9743
7 8
0.9891 0.9931
solution at t = 0.16
1 2 3 4 5 6
0.1000 0.3145 0.5094 0.6702 0.7905 0.8709
7 8
0.9159 0.9303
solution at t = 0.25
1 2 3 4 5 6
0.1000 0.2571 0.4052 0.5361 0.6434 0.7228
7 8
0.7713 0.7876
solution at t = 0.36
1 2 3 4 5 6
0.1000 0.2178 0.3295 0.4296 0.5129 0.5754
7 8
0.6142 0.6273
solution at t = 0.49
1 2 3 4 5 6
0.1000 0.1852 0.2661 0.3387 0.3993 0.4449
7 8
0.4732 0.4827
solution at t = 0.64
1 2 3 4 5 6
0.1000 0.1587 0.2144 0.2644 0.3061 0.3375
7 8
0.3570 0.3636
solution at t = 0.81
1 2 3 4 5 6
0.1000 0.1385 0.1752 0.2080 0.2354 0.2561
7 8
0.2689 0.2732
solution at t = 1.00
1 2 3 4 5 6
0.1000 0.1243 0.1475 0.1682 0.1855 0.1985
7 8
0.2066 0.2093
Example 2
Here, Problem C is solved from Sincovec and Madsen (1975). The equation is of diffusion-convection type with discontinuous coefficients. This problem illustrates a simple method for programming the evaluation routine for the derivative, ut. Note that the weak discontinuities at x = 0.5 are not evaluated in the expression for ut. The problem is defined as
#include <imsl.h>
#include <math.h>
#include <stdio.h>
void fcnut(int, float, float, float *, float *, float *, float *);
void fcnbc(int, float, float, float *, float *, float *);
int main()
{
int i, j, npdes = 1, nstep = 10, nx = 100;
float hinit, t = 0.0, tend, tol, xbreak[100], y[100];
char title[50];
void *state;
/* Set breakpoints and initial conditions */
for (i = 0; i < nx; i++) {
xbreak[i] = (float) i / (float) (nx - 1);
y[i] = 0.0;
}
y[0] = 1.0;
/* Initialize the solver */
tol = sqrt(imsl_f_machine(4));
hinit = 0.01*tol;
imsl_f_modified_method_of_lines_mgr(IMSL_PDE_INITIALIZE, &state,
IMSL_TOL, tol,
IMSL_HINIT, hinit,
0);
for (j = 1; j <= nstep; j++) {
tend = (float) j / (float) nstep;
/* Solve the problem */
imsl_f_modified_method_of_lines(npdes, &t, tend, nx, xbreak, y,
state, fcnut, fcnbc);
}
/* Print results at t=tend */
sprintf(title, "solution at t = %4.2f", t);
imsl_f_write_matrix(title, npdes, nx, y, 0);
}
void fcnut(int npdes, float x, float t, float *u, float *ux, float *uxx,
float *ut)
{
/* Define the PDE */
float d = 1.0, v = 1.0;
if (x <= 0.5) {
d = 5.0;
v = 1000.0;
}
ut[0] = d*uxx[0] - v*ux[0];
}
void fcnbc(int npdes, float x, float t, float *alpha, float *beta,
float *gam)
{
if (x == 0.0) {
*alpha = 1.0;
*beta = 0.0;
*gam = 1.0;
} else {
*alpha = 1.0;
*beta = 0.0;
*gam = 0.0;
}
}
Output
solution at t = 1.00
1 2 3 4 5 6
1.000 1.000 1.000 1.000 1.000 1.000
7 8 9 10 11 12
1.000 1.000 1.000 1.000 1.000 1.000
13 14 15 16 17 18
1.000 1.000 1.000 1.000 1.000 1.000
19 20 21 22 23 24
1.000 1.000 1.000 1.000 1.000 1.000
25 26 27 28 29 30
1.000 1.000 1.000 1.000 1.000 1.000
31 32 33 34 35 36
1.000 1.000 1.000 1.000 1.000 1.000
37 38 39 40 41 42
1.000 1.000 1.000 1.000 1.000 1.000
43 44 45 46 47 48
1.000 1.000 1.000 1.000 1.000 1.000
49 50 51 52 53 54
1.000 0.997 0.984 0.969 0.953 0.937
55 56 57 58 59 60
0.921 0.905 0.888 0.872 0.855 0.838
61 62 63 64 65 66
0.821 0.804 0.786 0.769 0.751 0.733
67 68 69 70 71 72
0.715 0.696 0.678 0.659 0.640 0.621
73 74 75 76 77 78
0.602 0.582 0.563 0.543 0.523 0.502
79 80 81 82 83 84
0.482 0.461 0.440 0.419 0.398 0.376
85 86 87 88 89 90
0.354 0.332 0.310 0.288 0.265 0.242
91 92 93 94 95 96
0.219 0.196 0.172 0.148 0.124 0.100
97 98 99 100
0.075 0.050 0.025 0.000
Example 3
In this example, using imsl_f_modified_method_of_lines, the linear normalized diffusion PDE ut = uxx is solved but with an optional use that provides values of the derivatives, ux, of the initial data. Due to errors in the numerical derivatives computed by spline interpolation, more precise derivative values are required when the initial data is u(x, 0) = 1 + cos[(2n − 1)πx], n > 1. The boundary conditions are “zero flux” conditions ux(0, t) = ux(1, t) = 0 for t > 0.
This optional usage signals that the derivative of the initial data is passed by the user. The values u(x, tend) and ux(x, tend) are output at the breakpoints with the optional usage.
#include <imsl.h>
#include <math.h>
#include <stdio.h>
void fcnut(int, float, float, float *, float *, float *, float *);
void fcnbc(int, float, float, float *, float *, float *);
int main()
{
int i, j, npdes = 1, nstep = 10, nx = 10;
float arg, hinit, tol, pi, t = 0.0, tend = 0.0;
float deriv[10], xbreak[10], y[10];
char title[50];
void *state;
pi = imsl_f_constant("pi", 0);
arg = 9.0 * pi;
/* Set breakpoints and initial conditions */
for (i = 0; i < nx; i++) {
xbreak[i] = (float) i / (float) (nx - 1);
y[i] = 1.0 + cos(arg * xbreak[i]);
deriv[i] = -arg * sin(arg * xbreak[i]);
}
/* Initialize the solver */
tol = sqrt(imsl_f_machine(4));
imsl_f_modified_method_of_lines_mgr(IMSL_PDE_INITIALIZE, &state,
IMSL_TOL, tol,
IMSL_INITIAL_VALUE_DERIVATIVE, deriv,
0);
for (j = 1; j <= nstep; j++) {
tend += 0.001;
/* Solve the problem */
imsl_f_modified_method_of_lines(npdes, &t, tend, nx, xbreak, y,
state, fcnut, fcnbc);
/* Print results at at every other t=tend */
if (!(j % 2)) {
sprintf(title, "\nsolution at t = %5.3f", t);
imsl_f_write_matrix(title, npdes, nx, y, 0);
sprintf(title, "\nderivative at t = %5.3f", t);
imsl_f_write_matrix(title, npdes, nx, deriv, 0);
}
}
}
void fcnut(int npdes, float x, float t, float *u, float *ux, float *uxx,
float *ut)
{
/* Define the PDE */
*ut = *uxx;
}
void fcnbc(int npdes, float x, float t, float *alpha, float *beta,
float *gam)
{
/* Define the boundary conditions */
alpha[0] = 0.0;
beta[0] = 1.0;
gam[0] = 0.0;
}
Output
solution at t = 0.002
1 2 3 4 5 6
1.234 0.766 1.234 0.766 1.234 0.766
7 8 9 10
1.234 0.766 1.234 0.766
derivative at t = 0.002
1 2 3 4 5 6
0.000e+000 8.983e-007 -3.682e-008 1.772e-006 4.368e-008 2.619e-006
7 8 9 10
-1.527e-006 4.956e-006 -3.003e-006 -3.009e-011
solution at t = 0.004
1 2 3 4 5 6
1.054 0.946 1.054 0.946 1.054 0.946
7 8 9 10
1.054 0.946 1.054 0.946
derivative at t = 0.004
1 2 3 4 5 6
0.000e+000 2.646e-007 4.763e-007 1.009e-006 -5.439e-007 -8.247e-007
7 8 9 10
3.142e-007 1.750e-006 -1.019e-006 4.300e-012
solution at t = 0.006
1 2 3 4 5 6
1.012 0.988 1.012 0.988 1.012 0.988
7 8 9 10
1.012 0.988 1.012 0.988
derivative at t = 0.006
1 2 3 4 5 6
0.000e+000 4.923e-007 4.082e-008 6.763e-007 -3.347e-007 -7.026e-007
7 8 9 10
4.525e-007 3.456e-007 -5.008e-007 1.327e-012
solution at t = 0.008
1 2 3 4 5 6
1.003 0.997 1.003 0.997 1.003 0.997
7 8 9 10
1.003 0.997 1.003 0.997
derivative at t = 0.008
1 2 3 4 5 6
0.000e+000 -1.323e-007 1.079e-006 2.271e-007 -7.651e-007 4.554e-007
7 8 9 10
7.479e-007 -5.015e-009 -3.918e-007 2.261e-013
solution at t = 0.010
1 2 3 4 5 6
1.001 0.999 1.001 0.999 1.001 0.999
7 8 9 10
1.001 0.999 1.001 0.999
derivative at t = 0.010
1 2 3 4 5 6
0.000e+000 9.523e-008 1.043e-006 3.912e-007 -6.791e-007 2.734e-008
7 8 9 10
4.506e-007 2.447e-007 -2.414e-008 -1.161e-015
Example 4
In this example, consider the linear normalized hyperbolic PDE, utt = uxx, the “vibrating string” equation. This naturally leads to a system of first order PDEs. Define a new dependent variable ut = v. Then, vt = uxx is the second equation in the system. Take as initial data u(x, 0) = sin(πx) and ut(x, 0) = v(x, 0) = 0. The ends of the string are fixed so u(0, t) = u(1, t) = v(0, t) = v(1, t) = 0. The exact solution to this problem is u(x, t) = sin(πx) cos(πt). Residuals are computed at the output values of t for 0 < t ≤ 2. Output is obtained at 200 steps in increments of 0.01.
Even though the sample code imsl_f_modified_method_of_lines gives satisfactory results for this PDE, users should be aware that for nonlinear problems, “shocks” can develop in the solution. The appearance of shocks may cause the code to fail in unpredictable ways. See Courant and Hilbert (1962), pp 488-490, for an introductory discussion of shocks in hyperbolic systems.
#include <imsl.h>
#include <math.h>
#include <stdio.h>
void fcnut(int, float, float, float *, float *, float *, float *);
void fcnbc(int, float, float, float *, float *, float *);
int main()
{
int i, j, npdes = 2, nstep = 200, nx = 10;
float error[10], erru = 0.0, err, hinit, pi, t = 0.0, tend = 0.0;
float deriv[20], tol, y[20], xbreak[10];
void *state;
pi = imsl_f_constant("pi", 0);
/* Set breakpoints and initial conditions */
for (i = 0; i < nx; i++) {
xbreak[i] = (float) i / (float) (nx - 1);
y[i] = sin(pi * xbreak[i]);
y[nx + i] = 0.0;
deriv[i] = pi * cos(pi * xbreak[i]);
deriv[nx + i] = 0.0;
}
/* Initialize the solver */
tol = sqrt(imsl_f_machine(4));
imsl_f_modified_method_of_lines_mgr(IMSL_PDE_INITIALIZE, &state,
IMSL_TOL, tol,
IMSL_INITIAL_VALUE_DERIVATIVE, deriv,
0);
for ( j=1; j <= nstep; j++) {
tend += 0.01;
/* Solve the problem */
imsl_f_modified_method_of_lines(npdes, &t, tend, nx, xbreak, y,
state, fcnut, fcnbc);
/* Look at output at steps of 0.01 and compute errors */
for (i = 0; i < nx; i++) {
error[i] = y[i] - sin(pi * xbreak[i]) * cos(pi *tend);
err = fabs(error[i]);
if (err > erru) erru = err;
}
}
printf("Maximum error in u(x,t) = %e\n", erru);
}
void fcnut(int npdes, float x, float t, float *u, float *ux, float *uxx,
float *ut)
{
/* Define the PDE */
ut[0] = u[1];
ut[1] = uxx[0];
}
void fcnbc(int npdes, float x, float t, float *alpha, float *beta,
float *gam)
{
/* Define the boundary conditions */
alpha[0] = 1.0;
beta[0] = 0.0;
gam[0] = 0.0;
alpha[1] = 1.0;
beta[1] = 0.0;
gam[1] = 0.0;
}
Output
Maximum error in u(x,t) = 5.757928e-003
Fatal Errors
|
IMSL_REPEATED_ERR_TEST_FAILURE |
After some initial success, the integration was halted by repeated error test failures. |
|
IMSL_INTEGRATION_HALTED_1
|
Integration halted after failing to pass the error test, even after reducing the stepsize by a factor of 1.0E+10 to H = #. TOL = # may be too small. |
|
IMSL_INTEGRATION_HALTED_2 |
Integration halted after failing to achieve corrector convergence, even after reducing the stepsize by a factor of 1.0E+10 to H = #. TOL = # may be too small. |
|
IMSL_INTEGRATION_HALTED_3 |
After some initial success, the integration was halted by a test on TOL = #. |
|
IMSL_TOL_TOO_SMALL_OR_STIFF |
On the next step X+H will equal X, with X = # and H = #. Either TOL = # is too small or the problem is stiff. |
