Chapter 5: Differential Equations

pde_method_of_lines

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.

Synopsis

#include <imsl.h>

void imsl_f_pde_method_of_lines_mgr (int task, void **state, ..., 0)

void imsl_f_pde_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_pde_method_of_lines_mgr and imsl_d_pde_method_of_lines.

Required Arguments for imsl_f_pde_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. These values for task are defined in the header 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_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 size 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_pde_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 gammap[])
User-supplied function to evaluate the boundary conditions. The boundary conditions accepted by imsl_f_pde_method_of_lines are

            Note: Users must supply the values ak and bk, which determine the values gk. Since gk can depend on t values of gk¢ also are required.

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 ak values.

float beta[]   (Output)
Array of length npdes containing the bk values.

float gammap[]   (Output)
Array of length npdes containing the derivatives,

Synopsis with Optional Arguments

#include <imsl.h>

void imsl_f_pde_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). 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 gammap[], 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 routine imsl_f_pde_method_of_lines uses the method of lines to solve the partial differential equation system

 

with the initial conditions

uk = uk(x, t)          at t = t0

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 fi(x) and yi(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, ai,k and bi,k. This system can be written in the matrixvector form as A dc/dt = F (t, y) 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 by differentiating the boundary conditions

for k = 1, ¼, M.

The initial conditions uk(x, t0) must satisfy the boundary conditions. Also, the
gk(t) must be continuous and have a smooth derivative, or the boundary conditions will not be properly imposed for t > t0.

If ak = bk = 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.

If the number of partial differential equations is M = 1 and the number of breakpoints is N = 4, then

The vector c is

c = [a1, b1, a2, b2, a3, b3, a4, b3]T

and the right-side F is

If M > 1, then each entry in the above matrix is replaced by an M ´ M diagonal matrix. The element a1 is replaced by diag(a1,1, ¼, a1,M). The elements aN, b1 and bN are handled in the same manner. The fi(pj) and yi(pj) elements are replaced by fi(pj)IM and yi(pj)IM where IM is the identity matrix of order M. See Madsen and Sincovec (1979) for further details about discretization errors and Jacobian matrix structure.

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 routine imsl_f_cub_spline_interp_e_cnd (Chapter 3, “Interpolation and Approximation”;;) to construct functions

such that

The IMSL routine imsl_f_cub_spline_value, Chapter 3, “Interpolation and Approximation;” is used to approximate the values

There is an optional use of imsl_f_pde_method_of_lines that 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. Some of the linear solvers were removed. Numerical Jacobians are used exclusively. The algorithm is unchanged. Gear's BDF method is used as the default because the system is typically stiff.

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 u1 = 0.1. This transition is completed by t = tδ = 0.09.

Due to restrictions in the type of boundary conditions successfully processed by imsl_f_pde_method_of_lines, it is necessary to provide the derivative boundary value function g¢ at x = 0 and at x = 1. The function g at x = 0 makes a smooth transition from the value u0 at t = t0 to the value u1 at t = tδ. The transition phase for g¢ is computed by evaluating a cubic interpolating polynomial. For this purpose, the function subprogram imsl_f_cub_spline_value, Chapter 3;, “Interpolation and Approximation” is used. The interpolation is performed as a first step in the user-supplied routine fcn_bc. The function and derivative values g(t0) = u0, g¢(t0) = 0, g(tδ) = u1, and g¢(tδ) = 0, are used as input to routine imsl_f_cub_spline_interp_e_cnd, to obtain the coefficients evaluated by imsl_f_cub_spline_value. Notice that g¢(t) = 0, t > tδ. The evaluation routine imsl_f_cub_spline_value will not yield this value so logic in the routine fcn_bc assigns g¢(t) = 0, t > tδ.

Link to example source

#include <imsl.h>

#include <math.h>

 

int main()

{

        void            fcnut(int, float, float, float *, float *, float *,

                                float *);

        void            fcnbc(int, float, float, float *, float *,

                                float *);

        int             npdes = 1;

        int             nx = 8;

        int             i;

        int             j = 1;

        int             nstep = 10;

        float           t = 0.0;

        float           tend;

        float           xbreak[8];

        float           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 */

 

        imsl_f_pde_method_of_lines_mgr(IMSL_PDE_INITIALIZE, &state,

                                       0);

 

        while (j <= nstep) {

                tend = (float) j++ / (float) nstep;

                tend *= tend;

 

                        /* Solve the problem */

 

                imsl_f_pde_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\0", 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 *gamp)

{

        static int      ndata;

        static int      first = 1;

        static float    delta = 0.09;

        static float    u0 = 1.0;

        static float    u1 = 0.1;

        static float    dfdata[2];

        static float    xdata[2];

        static float    fdata[2];

        static Imsl_f_ppoly *ppoly;

 

                        /* Compute interpolant first time only */

 

        if (first) {

                first = 0;

                ndata = 2;

                xdata[0] = 0.0;

                xdata[1] = delta;

                fdata[0] = u0;

                fdata[1] = u1;

                dfdata[0] = dfdata[1] = 0.0;

                ppoly = imsl_f_cub_spline_interp_e_cnd(ndata, xdata, fdata,

                                IMSL_LEFT, 1, dfdata[0],

                                IMSL_RIGHT, 1, dfdata[1],

                        0);

        }

 

                        /* Define boundary conditions */

 

        if (x == 0.0) {

 

                        /* These are for x = 0 */

 

                *alpha = 1.0;

                *beta = 0.0;

                *gamp = 0.0;

 

                        /* If in the boundary layer, compute

                           nonzero gamma prime */

 

                if (t <= delta)

                        *gamp = imsl_f_cub_spline_value(t, ppoly,

                                        IMSL_DERIV, 1,

                                        0);

        } else {

                        /* These are for x = 1 */

 

                *alpha = 0.0;

                *beta = 1.0;

                *gamp = 0.0;

        }

}

Output

                         solution at t = 0.01

         1           2           3           4           5           6

     0.969       0.997       1.000       1.000       1.000       1.000

 

         7           8

     1.000       1.000

 

                         solution at t = 0.04

         1           2           3           4           5           6

     0.625       0.871       0.962       0.991       0.998       1.000

          7           8

     1.000       1.000

 

                         solution at t = 0.09

         1           2           3           4           5           6

    0.1000      0.4602      0.7169      0.8671      0.9436      0.9781

 

         7           8

    0.9917      0.9951

 

                         solution at t = 0.16

         1           2           3           4           5           6

    0.1000      0.3130      0.5071      0.6681      0.7893      0.8708

 

         7           8

    0.9168      0.9315

 

                         solution at t = 0.25

         1           2           3           4           5           6

    0.1000      0.2567      0.4045      0.5354      0.6428      0.7224

 

         7           8

    0.7710      0.7874

 

                         solution at t = 0.36

         1           2           3           4           5           6

    0.1000      0.2176      0.3292      0.4292      0.5125      0.5751

 

         7           8

    0.6139      0.6270

 

                         solution at t = 0.49

         1           2           3           4           5           6

    0.1000      0.1852      0.2661      0.3386      0.3992      0.4448

 

         7           8

    0.4731      0.4827

 

                         solution at t = 0.64

         1           2           3           4           5           6

    0.1000      0.1588      0.2147      0.2648      0.3066      0.3381

 

         7           8

    0.3577      0.3643

 

                         solution at t = 0.81

         1           2           3           4           5           6

    0.1000      0.1387      0.1754      0.2083      0.2358      0.2565

 

         7           8

    0.2694      0.2738

 

                         solution at t = 1.00

         1           2           3           4           5           6

    0.1000      0.1242      0.1472      0.1678      0.1850      0.1980

 

         7           8

    0.2060      0.2087

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

Link to example source

#include <imsl.h>

#include <math.h>

 

int main()

{

        void            fcnut(int, float, float, float *, float *, float *,

                                float *);

        void            fcnbc(int, float, float, float *, float *,

                                float *);

        int             npdes = 1;

        int             nx = 100;

        int             i;

        int             j = 1;

        int             nstep = 10;

        float           t = 0.0;

        float           tend;

        float           xbreak[100];

        float           y[100];

        float           tol, hinit;

        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_pde_method_of_lines_mgr(IMSL_PDE_INITIALIZE, &state,

                                       IMSL_TOL, tol,

                                       IMSL_HINIT, hinit,

                                       0);

 

        while (j <= nstep) {

                tend = (float) j++ / (float) nstep;

 

                        /* Solve the problem */

 

                imsl_f_pde_method_of_lines(npdes, &t, tend, nx, xbreak, y,

                                           state, fcnut, fcnbc);

        }

                        /* Print results at t=tend */

 

                sprintf(title, "solution at t = %4.2f\0", 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 v;

        float d;

 

 

        if (x <= 0.5) {

                d = 5.0;

                v = 1000.0;

        }

        else

                d = v = 1.0;

 

        ut[0] = d*uxx[0] - v*ux[0];

}

 

void fcnbc(int npdes, float x, float t, float *alpha, float *beta,

        float *gamp)

{

        *alpha = 1.0;

        *beta = 0.0;

        *gamp = 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_pde_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)px], n > 1. The boundary conditions are “zero flux” conditions ux(0, t) = ux(1, t) = 0 for t > 0. Note that the initial data is compatible with these end conditions since the derivative function

vanishes at x = 0 and x = 1.

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.

Link to example source

#include <imsl.h>

#include <math.h>

 

int main()

{

        void            fcnut(int, float, float, float *, float *, float *,

                              float *);

        void            fcnbc(int, float, float, float *, float *, float *);

        int             npdes = 1;

        int             nx = 10;

        int             i;

        int             j = 1;

        int             nstep = 10;

        float           t = 0.0;

        float           tend = 0.0;

        float           xbreak[10];

        float           y[10], deriv[10];

        float           tol, hinit;

        float           pi, arg;

        char            title1[50];

        char            title2[50];

        void           *state;

 

        pi = imsl_d_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_pde_method_of_lines_mgr(IMSL_PDE_INITIALIZE, &state,

                                       IMSL_TOL, tol,

                                       IMSL_INITIAL_VALUE_DERIVATIVE,

                                       deriv,

                                       0);

 

        while (j <= nstep) {

                j++;

                tend += 0.001;

 

                        /* Solve the problem */

 

                imsl_f_pde_method_of_lines(npdes, &t, tend, nx, xbreak, y,

                                           state, fcnut, fcnbc);

 

                        /* Print results at at every other t=tend */

 

                if (j % 2) {

                        sprintf(title1, "\nsolution at t = %5.3f\0", t);

                        sprintf(title2, "\nderivative at t = %5.3f\0", t);

                        imsl_f_write_matrix(title1, npdes, nx, y, 0);

                        imsl_f_write_matrix(title2, npdes, nx, deriv, 0);

                }

        }

 

}

 

void fcnut(int npdes, float x, float t, float *u, float *ux, float *uxx,

      float *ut)

{

                        /* Define the PDE */

 

        ut[0] = uxx[0];

}

 

void fcnbc(int npdes, float x, float t, float *alpha, float *beta,

      float *gamp)

{

                        /* Define the boundary conditions */

 

        alpha[0] = 0.0;

        beta[0] = 1.0;

        gamp[0] = 0.0;

}

 Output

                         solution at t = 0.002

         1           2           3           4           5           6

     1.233       0.767       1.233       0.767       1.233       0.767

 

         7           8           9          10

     1.233       0.767       1.233       0.767

 

 

                        derivative at t = 0.002

         1           2           3           4           5           6

 0.000e+00  -5.172e-07   1.911e-06   1.818e-06  -5.230e-07   2.408e-06

 

         7           8           9          10

-2.517e-06   3.194e-06  -3.608e-06   2.023e-06

 

 

                         solution at t = 0.004

         1           2           3           4           5           6

     1.053       0.947       1.053       0.947       1.053       0.947

 

         7           8           9          10

     1.053       0.947       1.053       0.947

 

 

                        derivative at t = 0.004

         1           2           3           4           5           6

 0.000e+00  -1.332e-06  -9.059e-06  -4.401e-06   5.006e-06  -2.134e-06

 

         7           8           9          10

-1.733e-06   4.625e-06   6.741e-07   2.023e-06

 

 

 

                         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+00  -1.408e-06  -1.018e-06  -6.572e-07  -8.213e-07  -1.151e-06

 

         7           8           9          10

 1.051e-06   1.257e-06  -2.920e-07   2.023e-06

 

 

                         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+00  -1.028e-06   4.270e-06   3.114e-06  -3.085e-06  -1.492e-06

 

         7           8           9          10

 2.126e-06  -1.280e-06  -1.541e-06   2.023e-06

 

 

                         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+00  -7.596e-07   2.819e-07   1.547e-07  -1.469e-06  -9.516e-07

 

         7           8           9          10

 2.889e-07   8.956e-08   5.992e-07   2.023e-06

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(px) 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(px) cos(pt). 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_pde_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.

Link to example source

#include <imsl.h>

#include <math.h>

 

int main()

{

      void            fcnut(int, float, float, float *, float *, float *,

                            float *);

      void            fcnbc(int, float, float, float *, float *, float *);

      int             npdes = 2;

      int             nx = 10;

      int             i;

      int             j = 1;

      int             nstep = 200;

      float           t = 0.0;

      float           tend = 0.0;

      float           xbreak[20];

      float           y[20], deriv[20];

      float           tol, hinit;

      float           pi;

      float           error[10], erru;

      void           *state;

 

      pi = imsl_d_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_pde_method_of_lines_mgr(IMSL_PDE_INITIALIZE, &state,

                                     IMSL_TOL, tol,

                                     IMSL_INITIAL_VALUE_DERIVATIVE,

                                     deriv,

                                     0);

 

      while (j <= nstep) {

              j++;

              tend += 0.01;

                      /* Solve the problem */

 

              imsl_f_pde_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);

                      erru = imsl_f_max(erru, fabs(error[i]));

                }

        }

      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 *gamp)

{

                     /* Define the boundary conditions */

 

        alpha[0] = 1.0;

        beta[0] = 0.0;

        gamp[0] = 0.0;

        alpha[1] = 1.0;

        beta[1] = 0.0;

        gamp[1] = 0.0;

}

Output

Maximum error in u(x,t) = 6.228203e-04


Visual Numerics, Inc.
Visual Numerics - Developers of IMSL and PV-WAVE
http://www.vni.com/
PHONE: 713.784.3131
FAX:713.781.9260