pde1dMg

Solves a system of one-dimensional time-dependent partial differential equations using a moving-grid interface.

Synopsis

pde1dMgMgr (task, state)

pde1dMg (t, tend u, xl, xr, state, pdeSystems, boundaryConditions)

The function pde1dMgMgr is used to initialize and reset the problem, and the function pde1dMg is the integrator. The descriptions of both functions are provided below.

Required Arguments for pde1dMgMgr

int task (Input)

This function must be called with task set to PDE_INITIALIZE to set up for solving a system and with task equal to 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 pde1dMg

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 \(t_0\) where the integration of \(u_t\) begins. It is set to the value tend on return.

float tend (Input)

Mathematical value of \(t\) where the integration of \(u_t\) 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(x_j(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 pde1dMgMgr’s initialConditions.

float xl (Input)

Lower grid boundary, \(x_L\) .

float xr (Input)

Upper grid boundary, \(x_R\) .

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 pde1dMgMgr. It cannot be directly manipulated.

void pdeSystems(t, x, npdes, ngrids, full_u, grid_u, dudx, c, q, r, 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 \(m = 0\) is assumed, but this can be changed to one of the choices, \(m = 1, 2\). Use the optional arguments cartCoordinates, cylCoordinates, sphCoordinates for the respective values \(m = 0, 1, 2\). Return the values in the arrays as indicated:

\[\begin{split}\begin{array}{l} u^j = \mathtt{grid\_u}[j] \\ U = \mathtt{full\_u} \\ \frac{{\partial}u^j}{{\partial}x} = u_x^j = \mathtt{dudx}[j] \\ \text{c}[l][k] = C_{j,k}\left(x, t, u, u_x\right) \\ \text{r}[j] = r_j\left(x, t, u, u_x\right) \\ \text{q}[j] = q_j\left(x, t, u, u_x\right) \\ j, k = 0, \ldots, \mathit{NPDE} - 1 \\ \end{array}\end{split}\]

If any of the functions cannot be evaluated, set ires=3. Otherwise, do not change the value of ires.

void boundaryConditions (t, beta, gamma, full_u, grid_u, dudx, npdes, grids, left, ires) (Input)

User-supplied function to supply the boundary conditions, as expressed in Equation 2.

\[\begin{split}\begin{array}{l} u^j = \mathtt{grid\_u}[j] \\ U = \mathtt{full\_u} \\ \frac{\partial u^j}{\partial x} = u_x^j = \mathtt{dudx}[j] \\ \mathit{beta}[j] = \beta_j\left(x,t,u,u_x\right) \\ \mathtt{gamma}[j] = \gamma_j\left(x,t,u,u_x\right) \\ j = 0, \ldots \mathit{NPDE} - 1 \\ \end{array}\end{split}\]

The value \(x \in \left\{ x_L, x_R \right\}\), and the flag left=1 for \(x = x_L\). The flag has the value left=0 for \(x = x_R\). If any of the functions cannot be evaluated, set ires=3. Otherwise, do not change the value of ires.

Optional Arguments

cartCoordinates, or cylCoordinates, or sphCoordinates

cartCoordinates specifies cartesian coordinates, where \(m = 0\) in Equation 2. cylCoordinates specifies cylindrical or polar coordinates, where \(m = 1\) in Equation 2. sphCoordinates specifies spherical coordinates, where \(m = 2\) in Equation 2.

Default: cartCoordinates

timeSmoothing, float (Input)

Resets the value of the parameter \(\tau \geq 0\), described above.

Default: \(\tau = 0\).

spatialSmoothing, float (Input)

Resets the value of the parameter \(\kappa \geq 0\), described above.

Default: \(\kappa = 2\).

monitorRegularizing, float (Input)

Resets the value of the parameter \(\alpha \geq 0\), described above.

Default: \(\alpha = 0.01\).

maxBdfOrder, int maxBdfOrder, (Input)

Resets the maximum order for the bdf formulas used in differentialAlgebraicEqs. The new value can be any integer between 1 and 5. Some problems benefit by making this change. The default value of maxBdfOrder was chosen because differentialAlgebraicEqs may cycle on its selection of order and step-size with maxBdfOrder higher than value 2.

Default: maxBdfOrder=2.

userFactorSolve, int fac(neq, iband, a), void sol(neq, iband, g, 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.

initialConditions, void initialConditions(npdes, ngrids, 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 \(U_i \left( x_j(\mathit{tend}), \mathit{tend} \right)\) 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.

Optional Arguments

relativeTolerance, float (Input)

This option resets the value of the relative accuracy parameter used in differentialAlgebraicEqs.

Default: relativeTolerance=1.0E-2 for single precision, relativeTolerance=1.0E-4 for double precision.

absoluteTolerance, float (Input)

This option resets the value of the absolute accuracy parameter used in differentialAlgebraicEqs.

Default: absoluteTolerance=1E-2 for single precision, absoluteTolerance=1E-4 for double precision.

Examples

Remarks on the Examples

Due to its importance and the complexity of its interface, function pde1dMg 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.

Example 1 - Electrodynamics Model

This example is from Blom and Zegeling (1994). The system is

\[\begin{split}\begin{array}{l} u_t = {\varepsilon}pu_{xx} - g(u - v) \\ v_t = pv_{xx} + g(u - v), \\ \text{where } g(z) = \exp({\eta}z/3) - \exp(-2{\eta}z/3) \\ 0 \leq x \leq 1, 0 \leq t \leq 4 \\ u_x = 0 \text{ and } v = 0 \text{ at } x = 0 \\ u = 1 \text{ and } v_x = 0 \text{ at } x = 1 \\ \varepsilon = 0.143, p = 0.1743, \eta = 17.19 \\ \end{array}\end{split}\]

We make the connection between the model problem statement and the example:

\[\begin{split}\begin{array}{l} C = I_2 \\ m = 0, R_1 = \varepsilon pu_x, R_2 = pv_x \\ Q_1 = g(u-v), Q_2 = -Q_1 \\ u = 1 \text{ and } v = 0 \text{ at } t=1 \\ \end{array}\end{split}\]

The boundary conditions are

\[\begin{split}\begin{array}{l} \beta_1 = 1, \beta_2 = 0, \gamma_1 = 0, \gamma_2 = v, \text{ at } x = x_L = 0 \\ \beta_1 = 0, \beta_2 = 1, \gamma_1 = u - 1, \gamma_2 = 0, \text{ at } x = x_R = 1 \\ \end{array}\end{split}\]

This is a non-linear problem with sharply changing conditions near \(t = 0\). The default settings of integration parameters allow the problem to be solved. The use of pde1dMg requires two subroutines provided by the user to describe the differential equations, and boundary conditions.

from numpy import *
from pyimsl.math.pde1dMg import pde1dMg
from pyimsl.math.pde1dMgMgr import pde1dMgMgr, PDE_INITIALIZE, PDE_RESET


def initial_conditions(npdes, ngrids, u):
    for i in range(0, ngrids):
        u[0, i] = 1.0
        u[1, i] = 0.0


def pde_systems(t, x, npdes, ngrids, full_u, grid_u, dudx,
                c, q, r, ires):
    def setC(i, j, npdes, value):
        c[i * npdes + j] = value
    eps = 0.143
    eta = 17.19
    pp = 0.1743
    setC(0, 0, npdes, 1.0)
    setC(0, 1, npdes, 0.0)
    setC(1, 0, npdes, 0.0)
    setC(1, 1, npdes, 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]


def boundary_conditions(t, beta, gamma, full_u, grid_u, dudx,
                        ngrids, npdes, left, 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


npde = 2
nframes = 5
n = 51
u = zeros((npde + 1, n), dtype='double')
t0 = [0.0]
tout = 0.0
delta_t = 10.0
tend = 4.0
npdes = npde
ngrids = n
xl = 0.0
xr = 1.0

file1 = open("pde_ex01.out", mode="w")
file1.write("   %d\t%d\t%d" % (npdes, ngrids, nframes))
file1.write("\t%f\t%f\t%f\t%f\n" % (xl, xr, t0[0], tend))

# Initialize u
initial_conditions(npdes, ngrids, u)

state = []
pde1dMgMgr(PDE_INITIALIZE, state)
tout = 1.0e-3
while (True):
    pde1dMg(t0, tout, u, xl, xr, state,
            pde_systems, boundary_conditions)

    file1.write("%f\n" % (tout))
    for i in range(0, (npdes + 1)):
        for j in range(0, ngrids):
            file1.write("%16.10f    " % (u[i][j]))
            if (((j + 1) % 4) == 0):
                file1.write("\n")
        file1.write("\n")

    t0[0] = tout
    tout = tout * delta_t
    if (tend < tout):
        tout = tend
    if (t0[0] >= tend):
        break

file1.close()
pde1dMgMgr(PDE_RESET, state)

Example 2 - Inviscid Flow on a Plate

This example is a first order system from Pennington and Berzins, (1994). The equations are

\[\begin{split}\begin{array}{l} u_t = -v_x \\ uu_t = -vu_x + w_{xx} \\ w = u_x, \text{ implying that } uu_t = -vu_x + u_{xx} \\ u(0,t) = v(0,t) = 0,u(\infty,t) = u\left(x_R,t\right) = 1,t \geq 0 \\ u(x,0) = 1,v(x,0) = 0, x \geq 0 \\ \end{array}\end{split}\]

Following elimination of w, there remain \(\mathit{NPDE} = 2\) differential equations. The variable \(t\) is not time, but a second space variable. The integration goes from \(t = 0\) to \(t = 5\) . It is necessary to truncate the variable \(x\) at a finite value, say \(x_{\max} = x_R = 25\). In terms of the integrator, the system is defined by letting \(m = 0\) and

\[\begin{split}C = \left\{ C_{jk} \right\} = \begin{bmatrix} 10 \\ u0 \end{bmatrix}, R = \begin{bmatrix} -v \\ u_x \end{bmatrix}, Q = \begin{bmatrix} 0 \\ vu_x \end{bmatrix}\end{split}\]

The boundary conditions are satisfied by

\[\begin{split}\begin{array}{l} \beta = 0, \gamma = \begin{bmatrix} u = \exp(-20t) \\ v \end{bmatrix}, \text{ at } x = x_L \\ \\ \beta = 0, \gamma = \begin{bmatrix} u - 1 \\ v_x \end{bmatrix}, \text{ at } x = x_R \end{array}\end{split}\]

We use \(N = 10 + 51 = 61\) grid points and output the solution at steps of \(\mathit{\Delta}t = 0.1\).

This is a non-linear boundary layer problem with sharply changing conditions near \(t = 0\). The problem statement was modified so that boundary conditions are continuous near \(t = 0\). Without this change the underlying integration software, differentialAlgebraicEqs, cannot solve the problem. The continuous blending function \(u - \exp(-20t)\) is arbitrary and artfully chosen. This is a mathematical change to the problem, required because of the stated discontinuity at \(t = 0\). 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 \(X_L = 0\), 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.

from numpy import *
from pyimsl.math.pde1dMg import pde1dMg
from pyimsl.math.pde1dMgMgr import pde1dMgMgr, PDE_INITIALIZE, PDE_RESET

n1 = 10
n2 = 51
file1 = open("pde_ex02.out", mode="w")


def initial_conditions(npdes, ngrids, u):
    global n1, n2

    def setU(i, j, ngrids, value):
        u[i * ngrids + j] = value
    xl = 0.0
    xr = 25.0
    for i in range(0, ngrids):
        setU(0, i, ngrids, 1.0)
        setU(1, i, ngrids, 0.0)
        setU(2, i, ngrids, 0.0)
    dx1 = xr / n2
    dx2 = dx1 / n1
    # Grid
    for i in range(1, n1 + 1):
        i_ = i - 1
        setU(2, i_, ngrids, (i - 1) * dx2)
    for i in range(n1 + 1, n + 1):
        i_ = i - 1
        setU(2, i_, ngrids, (i - n1) * dx1)
    for i in range(0, npdes + 1):
        for j in range(0, ngrids):
            file1.write("%16.10f    " % (u[i * ngrids + j]))
            if (((j + 1) % 4) == 0):
                file1.write("\n")
        file1.write("\n")


def pde_systems(t, x, npdes, ngrids, full_u, grid_u,
                dudx, c, q, r, ires):
    def setC(i, j, npdes, value):
        c[i * npdes + j] = value
    setC(0, 0, npdes, 1.0)
    setC(1, 0, npdes, 0.0)
    setC(0, 1, npdes, grid_u[0])
    setC(1, 1, npdes, 0.0)
    r[0] = -grid_u[1]
    r[1] = dudx[0]
    q[0] = 0.0
    q[1] = grid_u[1] * dudx[0]


def boundary_conditions(t, beta, gamma, full_u, grid_u,
                        dudx, npdes, ngrids, left, ires):
    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]


npde = 2
n = n1 + n2
u = empty((npde + 1, n), dtype=double)
t0 = 0.0
delta_t = 1e-1
tend = 5.0
npdes = npde
ngrids = n
xl = 0.0
xr = 25.0
tau = 1.0e-3
atol = 1e-2
rtol = 0.0

nframes = int((tend + delta_t) / delta_t)
file1.write("   %d\t%d\t%d" % (npdes, ngrids, nframes))
file1.write("\t%f\t%f\t%f\t%f\n" % (xl, xr, t0, tend))
state = []

pde1dMgMgr(PDE_INITIALIZE, state,
           timeSmoothing=tau,
           initialConditions=initial_conditions)

t0 = [0.0]
tout = delta_t
while (True):

    pde1dMg(t0, tout, u, xl, xr, state,
            pde_systems, boundary_conditions,
            relativeTolerance=rtol,
            absoluteTolerance=atol)

    t0[0] = tout
    file1.write("%f\n" % (tout))
    for i in range(0, npdes + 1):
        for j in range(0, ngrids):
            file1.write("%16.10f    " % u[i][j])
            if (((j + 1) % 4) == 0):
                file1.write("\n")
        file1.write("\n")

    tout = tout + delta_t
    if (tend < tout):
        tout = tend
    if (t0[0] >= tend):
        break

file1.close()
pde1dMgMgr(PDE_RESET, state)

Example 3 - Population Dynamics

This example is from Pennington and Berzins (1994). The system is

\[\begin{split}\begin{array}{c} u_t = -u_x - I(t)u,x_L = 0 \leq x \leq a = x_R,t \geq 0 \hfill \\ I(t) = \int_0^a{u(x,t)dx} \hfill \\ u(x,0) = \frac{\exp(-x)}{2-\exp(-a)} \hfill \\ u(0,t) = g\left(\int_0^ab(x,I(t))u(x,t)dx,t\right), \text{where} \hfill \\ b(x,y) = \frac{xy\exp(-x)}{(y+1)^2}, \text{and} \hfill \\ g(z,t) = \hfill \\ \frac{4z(2-2\exp(-a)+\exp(-t))^2}{(1-\exp(-a))(1-(1+2a)\exp(-2a)(1-\exp(-a))+\exp(-t))} \end{array}\end{split}\]

This is a notable problem because it involves the unknown

\[u(x,t) = \frac{\exp(-x)}{1 - \exp(-a) + \exp(-t)}\]

across the entire domain. The software can solve the problem by introducing two dependent algebraic equations:

\[\begin{split}\begin{array}{l} v_1(t) = \int_0^a u(x,t) dx, \\ \\ v_2(t) = \int_0^a x \exp(-x) u(x,t) dx \\ \end{array}\end{split}\]

This leads to the modified system

\[\begin{split}\begin{array}{l} u_t = -u_x - v_1 u, \phantom{...} 0 \leq x \leq a, t \geq 0 \\ \\ u(0,t) = \tfrac{g(1,t)v_1v_2}{\left(v_1 + 1\right)^2} \\ \end{array}\end{split}\]

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 \(u(x,t)\) on the grid. The integrals are approximated using the trapezoid rule, commensurate with the truncation error in the integrator.

from numpy import *
from pyimsl.math.pde1dMg import pde1dMg
from pyimsl.math.pde1dMgMgr import pde1dMgMgr, PDE_INITIALIZE, PDE_RESET

file1 = open("pde_ex03.out", mode="w")


def initial_conditions(npdes, ngrids, u):
    def setU(i, j, ngrids, value):
        u[i * ngrids + j] = value
    global file1
    xl = 0.0
    xr = 5.0
    dx = (xr - xl) / (ngrids - 1)
    for i in range(0, ngrids):
        setU(0, i, ngrids, exp(-u[1 * ngrids + i]) / (2.0 - exp(-xr)))

    for i in range(0, npdes + 1):
        for j in range(0, ngrids):
            file1.write("%16.10f    " % (u[i * ngrids + j]))
            if (((j + 1) % 4) == 0):
                file1.write("\n")
        file1.write("\n")


def pde_systems(t, x, npdes, ngrids, full_u, grid_u,
                dudx, c, q, r, ires):
    sum = 0.0
    c[0] = 1.0
    r[0] = -1 * grid_u[0]
    for i in range(0, ngrids - 1):
        sum += (full_u[0 * ngrids + i] + full_u[0 * ngrids + (i + 1)]) * \
               (full_u[1 * ngrids + (i + 1)] - full_u[1 * ngrids + i])
    v1 = 0.5 * sum
    q[0] = v1 * grid_u[0]


def fcn_g(z, t):
    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


def boundary_conditions(t, beta, gamma, full_u, grid_u, dudx,
                        npdes, ngrids, left, ires):
    sum = 0.0
    sum1 = 0.0
    sum2 = 0.0
    sum3 = 0.0
    sum4 = 0.0

    for i in range(0, ngrids - 1):
        sum += (full_u[0 * ngrids + i]
                + full_u[0 * ngrids + (i + 1)]) * \
               (full_u[1 * ngrids + (i + 1)]
                - full_u[1 * ngrids + i])
        mid = 0.5 * (full_u[1 * ngrids + i] + full_u[1 * ngrids + (i + 1)])
        sum1 += mid * exp(-mid) * \
            ((full_u[0 * ngrids + i] + full_u[0 * ngrids + (i + 1)])
             * (full_u[1 * ngrids + (i + 1)] - full_u[1 * ngrids + 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]


npde = 1
n = 101
u = empty((npde + 1, n), dtype=double)
mid = array((n - 1), dtype=double)
npdes = npde
ngrids = n
t0 = [0.0]
tout = 0.0
delta_t = 1e-1
tend = 5.0
a = 5.0
state = []
xl = 0.0
xr = 5.0
tau = 1.0
atol = 1e-2
rtol = 0.0

nframes = (int)(tend + delta_t) / delta_t
file1.write("   %d\t%d\t%d" % (npdes, ngrids, nframes))
file1.write("\t%f\t%f\t%f\t%f\n" % (xl, xr, t0[0], tend))

pde1dMgMgr(PDE_INITIALIZE, state,
           timeSmoothing=tau,
           initialConditions=initial_conditions)

tout = delta_t
file1.write("%f\n" % (t0[0]))

while (True):
    pde1dMg(t0, tout, u, xl, xr, state,
            pde_systems, boundary_conditions,
            relativeTolerance=rtol,
            absoluteTolerance=atol)
    t0[0] = tout
    if (t0[0] <= tend):
        file1.write("%f\n" % (tout))
        for i in range(0, npdes + 1):
            for j in range(0, ngrids):
                file1.write("%16.10f    " % (u[i][j]))
                if (((j + 1) % 4) == 0):
                    file1.write("\n")
            file1.write("\n")
    if ((tout + delta_t) < tend):
        tout = tout + delta_t
    else:
        tout = tend

    if (t0[0] >= tend):
        break

file1.close()
pde1dMgMgr(PDE_RESET, state)

Example 4 - A Model in Cylindrical Coordinates

This example is from Blom and Zegeling (1994). The system models a reactor-diffusion problem:

\[\begin{split}\begin{array}{l} T_z = r^{-1}\frac{\partial\left({\beta}rT_r\right)}{{\partial}r} + \gamma\exp\left(\frac{T}{1+{\varepsilon}T}\right) \\ T_r(0,z) = 0, T(1,z) = 0, z > 0 \\ T(r,0) = 0, 0 \leq r < 1 \\ \beta = 10^{-4}, \gamma = 1, \varepsilon = 0.1 \\ \end{array}\end{split}\]

The axial direction \(z\) is treated as a time coordinate. The radius \(r\) is treated as the single space variable.

This is a non-linear problem in cylindrical coordinates. Our example illustrates assigning \(m = 1\) in Equation 2. We provide the optional argument cylCoordinates that resets this value from its default, \(m = 0\).

from numpy import *
from pyimsl.math.pde1dMg import pde1dMg
from pyimsl.math.pde1dMgMgr import pde1dMgMgr, PDE_INITIALIZE, PDE_RESET


def initial_conditions(npdes, ngrids, t):
    for i in range(0, ngrids):
        t[0][i] = 0.0


def pde_systems(t, x, npdes, ngrids, u, grid_u, dudx,
                c, q, r, ires):
    beta = 01e-4
    gamma = 1.0
    eps = 1e-1

    c[0] = 1.0
    r[0] = beta * dudx[0]
    q[0] = -1.0 * gamma * exp(grid_u[0] / (1.0 + eps * grid_u[0]))


def boundary_conditions(t, beta, gamma, u, grid_u, dudx,
                        npdes, ngrids, left, ires):
    if (left):
        beta[0] = 1.0
        gamma[0] = 0.0
    else:
        beta[0] = 0.0
        gamma[0] = grid_u[0]


npde = 1
n = 41
t = empty((npde + 1, n), dtype=double)
z0 = [0.0]
delta_z = 1e-1
zend = 1.0
zmax = 1.0
beta = 1e-4
gamma = 1.0
eps = 1e-1
npdes = npde
ngrids = n
xl = 0.0
xr = 1.0
m = 1

file1 = open("pde_ex04.out", mode="w")
nframes = (int)((zend + delta_z) / delta_z) - 1
file1.write("   %d\t%d\t%d" % (npdes, n, nframes))
file1.write("\t%f\t%f\t%f\t%f\n" % (xl, xr, z0[0], zend))

state = []
pde1dMgMgr(PDE_INITIALIZE, state, cylCoordinates=True)
initial_conditions(npdes, ngrids, t)

zout = delta_z
while (True):
    pde1dMg(z0, zout, t, xl, xr, state,
            pde_systems, boundary_conditions)

    z0[0] = zout
    if (z0[0] <= zend):
        file1.write("%f\n" % (zout))
        for i in range(0, npdes + 1):
            for j in range(0, ngrids):
                file1.write("%16.10f    " % (t[i, j]))
                if (((j + 1) % 4) == 0):
                    file1.write("\n")
            file1.write("\n")

    if ((zout + delta_z) < zend):
        zout = zout + delta_z
    else:
        zout = zend

    if (z0[0] >= zend):
        break

file1.close()
pde1dMgMgr(PDE_RESET, state)

Example 5 - A ‘Hot Spot’ Model

This example is presented more fully in Verwer, et al., (1989). The system is a normalized problem relating the temperature \(u(x,t)\), of a reactant in a chemical system. The formula for \(h(z)\) is equivalent to their example.

\[\begin{split}\begin{array}{l} u_t = u_{xx} + h(u), \\ \text{where } h(z) = \tfrac{R}{a\delta}(1 + a - z)\exp(-\delta(1/z-1)), \\ a = 1,\delta = 20, R = 5 \\ 0 \leq x \leq 1,0 \leq t \leq 0.29 \\ u(x,0) = 1 \\ u_x = 0,x = 0 \\ u = 1,x = 1 \\ \end{array}\end{split}\]

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.

from numpy import *
from pyimsl.math.pde1dMg import pde1dMg
from pyimsl.math.pde1dMgMgr import pde1dMgMgr, PDE_INITIALIZE, PDE_RESET


def initial_conditions(npdes, ngrids, u):
    for i in range(0, ngrids):
        u[0][i] = 1.0


def pde_systems(t, x, npdes, ngrids, full_u, grid_u, dudx,
                c, q, r, ires):
    c[0] = 1.0
    r[0] = dudx[0]
    q[0] = -fcn_h(grid_u[0])


def boundary_conditions(t, beta, gamma, full_u, grid_u, dudx,
                        npdes, ngrids, left, ires):
    if (left):
        beta[0] = 0.0
        gamma[0] = dudx[0]
    else:
        beta[0] = 0.0
        gamma[0] = grid_u[0] - 1.0


def fcn_h(z):
    a = 1.0
    delta = 2e1
    r = 5.0

    return (r / (a * delta)) * (1.0 + a - z) * \
        exp(-delta * (1.0 / z - 1.0))


npde = 1
n = 80
u = empty((npde + 1, n), dtype=double)
t0 = [0.0]
tout = 0.0
delta_t = 1e-2
tend = 29e-2
u0 = 1.0
u1 = 0.0
tdelta = 1e-1
tol = 29e-2
a = 1.0
delta = 20.0
r = 5.0
state = []
npdes = npde
ngrids = n
xl = 0.0
xr = 1.0
max_bdf_order = 5

file1 = open("pde_ex05.out", mode="w")
nframes = int((tend + delta_t) / delta_t) - 1
file1.write("   %d\t%d\t%d" % (npdes, ngrids, nframes))
file1.write("\t%f\t%f\t%f\t%f\n" % (xl, xr, t0[0], tend))

initial_conditions(npdes, ngrids, u)

pde1dMgMgr(PDE_INITIALIZE, state, maxBdfOrder=max_bdf_order)

tout = delta_t

while (True):
    pde1dMg(t0, tout, u, xl, xr, state,
            pde_systems, boundary_conditions)

    t0[0] = tout
    if (t0[0] <= tend):
        file1.write("%f\n" % (tout))
        for i in range(0, npdes + 1):
            for j in range(0, ngrids):
                file1.write("%16.10f    " % (u[i][j]))
                if (((j + 1) % 4) == 0):
                    file1.write("\n")
    if ((tout + delta_t) < tend):
        tout = tout + delta_t
    else:
        tout = tend

    if (t0[0] >= tend):
        break

file1.close()
pde1dMgMgr(PDE_RESET, state)

Example 6 - Traveling Waves

This example is presented more fully in Verwer, et al., (1989). The system is a normalized problem relating the interaction of two waves, \(u(x,t)\) and \(v(x,t)\) 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.

\[\begin{split}\begin{array}{c} u_t = -u_x-100uv, \hfill \\ v_t = v_x-100uv, \hfill \\ -0.5 \leq x \leq 0.5,0 \leq t \leq 0.5 \hfill \\ \begin{aligned} u(x,0) &= 0.5(1+\cos(10{\pi}x)), x \in \left[-0.3,-0.1\right], \text{and} \\ &= 0, \text{otherwise}, \end{aligned} \hfill \\ \begin{aligned} v(x,0) &= 0.5(1+\cos(10{\pi}x)), x \in \left[0.1,0.3\right], \text{and} \\ &=0, \text{otherwise}, \end{aligned} \hfill \\ u = v = 0 \text{ at both ends}, t \geq 0 \hfill \\ \end{array}\end{split}\]

This is a non-linear system of first order equations.

from numpy import *
from pyimsl.math.pde1dMg import pde1dMg
from pyimsl.math.pde1dMgMgr import pde1dMgMgr, PDE_INITIALIZE, PDE_RESET

file1 = open("pde_ex06.out", mode="w")


def initial_conditions(npdes, ngrids, u):
    def setU(ngrids_, i_, j_, value_):
        u[i_ * ngrids_ + j_] = value_

    def getU(ngrids_, i_, j_):
        return u[i_ * ngrids_ + j_]

    for i in range(0, ngrids):
        pulse = (0.5 * (1.0 + cos(10.0 * pi * getU(ngrids, npdes, i))))

        setU(ngrids, 0, i, pulse)
        setU(ngrids, 1, i, pulse)

    for i in range(0, ngrids):
        utmp = getU(ngrids, npdes, i)
        if ((utmp < -3e-1) or (utmp > -1e-1)):
            setU(ngrids, 0, i, 0.0)
        if ((utmp < 1e-1) or (utmp > 3e-1)):
            setU(ngrids, 1, i, 0.0)
    for i in range(0, npdes + 1):
        for j in range(0, ngrids):
            utmp = getU(ngrids, i, j)
            file1.write("%16.10f    " % (utmp))
            if (((j + 1) % 4) == 0):
                file1.write("\n")
        file1.write("\n")


def pde_systems(t, x, npdes, ngrids, full_u, grid_u, dudx,
                c, q, r, ires):
    def setC(npdes_, i_, j_, value_):
        c[i_ * npdes_ + j_] = value_
    setC(npdes, 0, 0, 1.0)
    setC(npdes, 0, 1, 0.0)
    setC(npdes, 1, 0, 0.0)
    setC(npdes, 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]


def boundary_conditions(t, beta, gamma, full_u, grid_u, dudx,
                        npdes, ngrids, left, ires):
    beta[0] = 0.0
    beta[1] = 0.0
    gamma[0] = grid_u[0]
    gamma[1] = grid_u[1]


npde = 2
n = 50
npdes = npde
ngrids = n
xl = -0.5
xr = 0.5

u = empty((npde + 1, n), dtype=double)
t0 = [0.0]
tout = 0.0
delta_t = 5e-2
tend = 5e-1
state = []
tau = 1e-3
atol = 1e-3
rtol = 0.0
max_bdf_order = 3

nframes = int((tend + delta_t) / delta_t)
file1.write("   %d\t%d\t%d" % (npdes, ngrids, nframes))
file1.write("\t%f\t%f\t%f\t%f\n" % (xl, xr, t0[0], tend))

pde1dMgMgr(PDE_INITIALIZE, state,
           timeSmoothing=tau,
           maxBdfOrder=max_bdf_order,
           initialConditions=initial_conditions)

file1.write("%f\n" % (t0[0]))
tout = delta_t

while (True):
    pde1dMg(t0, tout, u, xl, xr, state,
            pde_systems, boundary_conditions,
            relativeTolerance=rtol,
            absoluteTolerance=atol)
    t0[0] = tout
    if (t0[0] <= tend):
        file1.write("%f\n" % (tout))
        for i in range(0, npdes + 1):
            for j in range(0, ngrids):
                file1.write("%16.10f    " % (u[i, j]))
                if (((j + 1) % 4) == 0):
                    file1.write("\n")
            file1.write("\n")

    if ((tout + delta_t) < tend):
        tout = tout + delta_t
    else:
        tout = tend

    if (t0[0] >= tend):
        break

file1.close()
pde1dMgMgr(PDE_RESET, state)

Example 7 - Black-Scholes

The value of a European “call option,” \(c(s,t)\), with exercise price \(e\) and expiration date \(T\), satisfies the “asset-or-nothing payoff ” \(c(s,T) = s, s \geq e; = 0, s < e\). Prior to expiration \(c(s,t)\) is estimated by the Black-Scholes differential equation \(c_t + \tfrac{\sigma^2}{2} s^2 c_{ss} + rsc_s - rc \equiv c_t + \tfrac{\sigma^2}{2} \left( s^2 c_s \right)_s + \left( r - \sigma^2 \right) sc_s - rc = 0\) . The parameters in the model are the risk-free interest rate, \(r\), and the stock volatility, \(\sigma\). The boundary conditions are \(c(0,t) = 0\) and \(c_s(s,t) \approx 1,s \rightarrow \infty\). 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 PyIMSL function normalCdf, see Chapter 9, “Special Functions.” With numerical integration the equation itself or the payoff can be readily changed to include other formulas, \(c(s, T)\) , and corresponding boundary conditions. We use \(e = 100, r = 0.08, T -t = 0.25, \sigma^2 = 0.04, s_L = 0 \text{ and } s_R = 150\) .

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 \(10^{-3}\). In $US, this is one-tenth of a cent.

from numpy import *
from pyimsl.math.pde1dMg import pde1dMg
from pyimsl.math.pde1dMgMgr import pde1dMgMgr, PDE_INITIALIZE, PDE_RESET


def initial_conditions(npdes, ngrids, u):
    xl = 0.0
    xr = 150.0
    e = 100.0
    dx = (xr - xl) / (ngrids - 1)
    for i in range(0, ngrids):
        xi = xl + i * dx
        if (xi <= e):
            u[0][i] = 0.0
        else:
            u[0][i] = xi


def pde_systems(t, x, npdes, ngrids, full_u, grid_u, dudx, c,
                q, r, ires):
    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]


def boundary_conditions(t, beta, gamma, full_u, grid_u, dudx,
                        npdes, ngrids, left, ires):
    if (left):
        beta[0] = 0.0
        gamma[0] = grid_u[0]
    else:
        beta[0] = 0.0
        gamma[0] = dudx[0] - 1.0


npde = 1
n = 100
xl = 0.0
xr = 150.0

u = empty((npde + 1, n), dtype=double)
t0 = [0.0]
delta_t = 25e-3
tend = 25e-2
xmax = 150.0
state = []
npdes = npde
ngrids = n
tau = 5e-3
atol = 1e-2
rtol = 0.0
max_bdf_order = 5

file1 = open("pde_ex07.out", mode="w")

nframes = int((tend + delta_t) / delta_t)
file1.write("   %d\t%d\t%d" % (npdes, ngrids, nframes))
file1.write("\t%f\t%f\t%f\t%f\n" % (xl, xr, t0[0], tend))

initial_conditions(npdes, ngrids, u)

pde1dMgMgr(PDE_INITIALIZE, state,
           timeSmoothing=tau,
           maxBdfOrder=max_bdf_order)

tout = delta_t

while (True):
    pde1dMg(t0, tout, u, xl, xr, state,
            pde_systems, boundary_conditions,
            relativeTolerance=rtol,
            absoluteTolerance=atol)

    t0[0] = tout
    if (t0[0] <= tend):
        file1.write("%f\n" % (tout))
        for i in range(0, npdes + 1):
            for j in range(0, ngrids):
                file1.write("%16.10f    " % (u[i][j]))
                if (((j + 1) % 4) == 0):
                    file1.write("\n")
    if ((tout + delta_t) < tend):
        tout = tout + delta_t
    else:
        tout = tend

    if (t0[0] >= tend):
        break

file1.close()
pde1dMgMgr(PDE_RESET, state)

Fatal Errors

IMSL_STOP_USER_FCN

Request from user supplied function to stop algorithm. User flag = “#”.