modifiedMethodOfLines¶

Solves a system of partial differential equations of the form \(u_t = f(x, t, u, u_x, u_{xx})\) using the method of lines. The solution is represented with cubic Hermite polynomials.

Synopsis¶

modifiedMethodOfLinesMgr (task, state)

modifiedMethodOfLines (npdes, t, tend, nx, xbreak, y, state, fcnUt, fcnBc)

Required Arguments for modifiedMethodOfLinesMgr¶

int task (Input): This function must be called with task set to PDE_INITIALIZE to set up memory and default values prior to solving a problem and with task equal to 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 modifiedMethodOfLines¶

int npdes (Input): Number of differential equations.
float t (Input/Output): Independent variable. On input, t supplies the initial time, \(t_0\). 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] = \(u_k\)(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 modifiedMethodOfLinesMgr. It cannot be directly manipulated.
void fcnUt(npdes, x, t, u[], ux[], uxx[], ut[]): User-supplied function to evaluate \(u_t\).
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, \(u_x\).
float uxx[] (Input): Array of length npdes containing the second derivative, \(u_{xx}\).
float ut[] (Output): Array of length npdes containing the computed derivatives \(u_t\).
void fcnBc(npdes, x, t, alpha[], beta[], gamma[]): User-supplied function to evaluate the boundary conditions. The boundary conditions accepted by modifiedMethodOfLines are

\[\alpha_k u_k + \beta_k \frac{\partial u_k}{\partial x} = \gamma_k\]

NOTE: Users must supply the values \(\alpha_k\), \(\beta_k\), and \(y_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 \(\alpha_k\) values.
float beta[] (Output): Array of length npdes containing the \(\beta_k\) values.
float gamma[] (Output): Array of length npdes containing the values of \(y_k\).

Optional Arguments¶

tol, float (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\*machine(4)

initialStepsize, float (Input)

Initial step size in the t integration. This value must be nonnegative. If initialStepsize is zero, an initial step size of 0.001|tend - \(t_0\)| will be arbitrarily used. The step will be applied in the direction of integration.

Default: initialStepsize = 0.0

initialValueDerivative, float [[]] (Input/Output)

Supply the derivative values \(u_x(x, t_0)\) in initialValueDerivative, an array of length npdes by nx. This derivative information is input as

\[\mathrm{initialValueDerivative[k,i]} = \frac{\partial u_k}{\partial x} \left(x, t_0\right) \text{ at } x = x[i]\]

The array initialValueDerivative contains the derivative values as output:

\[\mathrm{initialValueDerivative[k,i]} = \frac{\partial u_k}{\partial x} \left(x, \mathrm{tend}\right) \text{ at } x = x[i]\]

Default: Derivatives are computed using cubic spline interpolation.

htrial (Output)

Return the current trial step size.

Description¶

Let M = npdes, N = nx and \(x_i\) = xbreak(I). The function modifiedMethodOfLines uses the method of lines to solve the partial differential equation system

\[\mathrm{initialValueDerivative[k,i]} = \frac{\partial u_k}{\partial x} \left(x, t_0\right) \text{ at } x = x[i]\]

\[\mathrm{initialValueDerivative[k,i]} = \frac{\partial u_k}{\partial x} \left(x, \mathrm{tend}\right) \text{ at } x = x[i]\]

\[\frac{\partial u_k}{\partial t} = f_k \left( x,t,u_1,\ldots u_M, \frac{\partial u_1}{\partial x}, \ldots \frac{\partial u_M}{\partial x}, \frac{\partial^2 u_1}{\partial x^2}, \ldots \frac{\partial^2 u_M}{\partial x^2} \right)\]

with the initial conditions

\[u_k = u_k(x,t) \text{ at } t = t_0\]

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

\[\alpha_k u_k + \beta_k \frac{\partial u_k}{\partial x} = \gamma_k \text{ at } x = x_1 \text{ and at } x = x_N\]

\[\hat{u}_k(x,t) = \sum_{i=1}^{N} \left(a_{k,i}(t) \phi_i(x) + b_{k,i}(t) \psi_i(x) \right)\]

where \(\phi_i(x)\) and \(\psi_i(x)\) are the standard basis functions for the cubic Hermite polynomials with the knots \(x_1 < x_2 < ... < x_N\). These are piecewise cubic polynomials with continuous first derivatives. At the breakpoints, they satisfy

\[\begin{split}\begin{array}{c} \phi_i\left(x_l\right) = \delta_{il} \phantom{...} \psi_i\left(x_l\right) = 0 \\ \tfrac{d\phi_i}{dx}\left(x_l\right) = 0 \phantom{...} \tfrac{d\psi_i}{dx}\left(x_l\right) = \delta_{il} \end{array}\end{split}\]

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,

\[\begin{split}\begin{array}{c} p_{2j-1} = x_j + \frac{3-\sqrt{3}}{6} \left(x_{j+1} - x_j\right) \\ p_{2j} = x_j + \frac{3 + \sqrt{3}}{6} \left(x_{j+1} - x_j\right) \\ \end{array}\end{split}\]

for \(j = 1, ..., N\). The collocation approximation to the differential equation is

\[\begin{split}\begin{array}{l} \displaystyle\sum_{i=1}^N \frac{da_{k,i}}{dt}\phi_i\left(p_j\right) + \frac{db_{k,i}}{dt}\psi_i\left(p_j\right) = \\ f_k\left(p_j,t,\hat{u}_1\left(p_j\right), \ldots, \hat{u}_M\left(p_j\right), \dots, \left(\hat{u}_1\right)_{xx}\left(p_j\right), \ldots, \left(\hat{u}_M\right)_{xx}\left(p_j\right)\right) \end{array}\end{split}\]

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, \(a_{k,i}\) and \(b_{k,i}\). This system can be written in the matrix‑vector form as \(A dc / dt = F(t, c)\) with \(c(t_0) = c_0\) where c is a vector of coefficients of length 2M N and \(c_0\) 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 \(a_{k,i}\). The initial values of the \(b_{k,i}\) are obtained by using the PyIMSL cubic spline function cubSplineInterpECnd (Interpolation and Approximation) to construct functions

\[\hat{u}_k \left( x_i, t_0 \right)\]

such that

\[d\hat{u}_k \left( x_i, t_0 \right) = a_{ki}\]

The PyIMSL function cubSplineValue (Interpolation and Approximation) is used to approximate the values

\[\hat{u}_k \left(x_i, t_0\right) = a_{ki}\]

\[\frac{d\hat{u}_k}{dx} \left(x_i, t_0\right) \equiv b_{k,i}\]

Optional argument initialValueDerivative allows the user to provide the initial values of \(b_{k,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 odeAdamsKrogh. 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 PyIMSL 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, \(u_t = u_{xx}\), \(0 \leq x \leq 1\), \(t > t_0\), is solved. The initial values are \(t_0 = 0\), \(u(x, t_0) = u_0 = 1\). There is a “zero-flux” boundary condition at \(x = 1\), namely \(u_x(1, t) = 0\), (\(t > t_0\)). The boundary value of \(u(0, t)\) is abruptly changed from \(u_0\) 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.

from numpy import *
from pyimsl.math.modifiedMethodOfLines import modifiedMethodOfLines
from pyimsl.math.modifiedMethodOfLinesMgr import modifiedMethodOfLinesMgr, PDE_INITIALIZE, PDE_RESET
from pyimsl.math.writeMatrix import writeMatrix


def fcnut(npdes, x, t, u, ux, uxx, ut):
    # Define the PDE
    ut[0] = uxx[0]


def fcnbc(npdes, x, t, alpha, beta, gam):
    delta = 0.09
    u0 = 1.0
    u1 = 0.1
    # Define boundary conditions
    if (x == 0.0):
        # These are for x=0
        alpha[0] = 1.0
        beta[0] = 0.0
        gam[0] = u1
        # If in the boundary layer, compute
        # nonzero gamma
        if (t <= delta):
            gam[0] = u0 + (u1 - u0) * t / delta
    else:
        # These are for x=1
        alpha[0] = 0.0
        beta[0] = 1.0
        gam[0] = 0.0


npdes = 1
nx = 8
nstep = 10
t = [0.0]
tend = 0.0
xbreak = zeros((8), dtype='double')
y = zeros((1, 8), dtype='double')
state = []

# Set breakpoints and initial conditions
for i in range(0, nx):
    xbreak[i] = double(i) / double(nx - 1)
    y[0, i] = 1.0

# Initialize the solver
tol = 10.e-4
hinit = 0.01 * tol
modifiedMethodOfLinesMgr(PDE_INITIALIZE, state,
                         tol=tol, hinit=hinit)
# while (j <= nstep):
for j in range(1, nstep + 1):
    tend = float(j) / float(nstep)
    tend *= tend

    # Solve the problem
    modifiedMethodOfLines(t, float(tend), xbreak, y, state, fcnut, fcnbc)

    # Print results at current t=tend
    title = "Solution at t = %4.2f" % (tend)
    writeMatrix(title, y)

modifiedMethodOfLinesMgr(PDE_RESET, state)

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.1002       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.3143       0.5092       0.6701       0.7906       0.8710
 
          7            8
     0.9160       0.9304
 
                            Solution at t = 0.25
          1            2            3            4            5            6
     0.1000       0.2571       0.4052       0.5361       0.6435       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.5755
 
          7            8
     0.6142       0.6273
 
                            Solution at t = 0.49
          1            2            3            4            5            6
     0.1000       0.1853       0.2662       0.3389       0.3995       0.4451
 
          7            8
     0.4734       0.4830
 
                            Solution at t = 0.64
          1            2            3            4            5            6
     0.1000       0.1589       0.2149       0.2650       0.3070       0.3385
 
          7            8
     0.3580       0.3647
 
                            Solution at t = 0.81
          1            2            3            4            5            6
     0.1000       0.1388       0.1756       0.2086       0.2362       0.2569
 
          7            8
     0.2698       0.2742
 
                            Solution at t = 1.00
          1            2            3            4            5            6
     0.1000       0.1243       0.1474       0.1680       0.1853       0.1983
 
          7            8
     0.2064       0.2091

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, \(u_t\). Note that the weak discontinuities at \(x = 0.5\) are not evaluated in the expression for \(u_t\). The problem is defined as

\[\begin{split}\begin{array}{c} u_t = \partial u/ \partial t = \partial / \partial x (D(x) \partial u / \partial x) - v(x) \partial u / \partial x \\ x \in [0,1],t > 0 \end{array}\end{split}\]

\[\begin{split}D(x) = \begin{cases} 5 & \text{if } 0 \leq x < 0.5 \\ 1 & \text{if } 0.5 < x \leq 1.0 \\ \end{cases}\end{split}\]

\[\begin{split}\begin{array}{c} \nu(x) = \begin{cases} 1000.0 & \text{if } 0 \leq x < 0.5 \\ 1 & \text{if } 0.5 < x \leq 1.0 \\ \end{cases} \\ u(x,0) = \begin{cases} 1 & \text{if } x = 0 \\ 0 & \text{if } x > 0 \\ \end{cases} \\ u(0,t) = 1,u(1,t) = 0 \\ \end{array}\end{split}\]

from numpy import *
from pyimsl.math.machine import machine
from pyimsl.math.modifiedMethodOfLines import modifiedMethodOfLines
from pyimsl.math.modifiedMethodOfLinesMgr import modifiedMethodOfLinesMgr, PDE_INITIALIZE, PDE_RESET
from pyimsl.math.writeMatrix import writeMatrix


def fcnut(npdes, x, t, u, ux, uxx, ut):
    # Define the PDE
    if (x <= 0.5):
        d = 5.0
        v = 1000.0
    else:
        d = v = 1.0
    ut[0] = d * uxx[0] - v * ux[0]


def fcnbc(npdes, x, t, alpha, beta, gam):
    alpha[0] = 1.0
    beta[0] = 0.0
    gam[0] = 0.0
    if (x == 0):
        gam[0] = 1.0


npdes = 1
nx = 100
j = 1
nstep = 10
t = [0.0]
tend = 0.0
xbreak = zeros((100), dtype='double')
y = zeros((1, 100), dtype='double')
state = []

# Set breakpoints and initial conditions
for i in range(0, nx):
    xbreak[i] = double(i) / double(nx - 1)
    y[0, i] = 0.0
y[0, 0] = 1.0

# Initialize the solver
tol = sqrt(machine(4))
hinit = 0.01 * tol
modifiedMethodOfLinesMgr(PDE_INITIALIZE, state,
                         tol=tol, hinit=hinit)
for j in range(1, nstep + 1):
    tend = float(j) / float(nstep)
    # Solve the problem
    modifiedMethodOfLines(t, tend, xbreak, y, state, fcnut, fcnbc)

# Print results at t=tend
title = "Solution at t = %4.2f" % (tend)
writeMatrix(title, y)
modifiedMethodOfLinesMgr(PDE_RESET, state)

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 modifiedMethodOfLines, the linear normalized diffusion PDE \(u_t = u_{xx}\) is solved but with an optional use that provides values of the derivatives, \(u_x\), 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)\pi x]\), \(n > 1\). The boundary conditions are “zero flux” conditions \(u_x(0, t) = u_x(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 \(u_x(x, tend)\) are output at the breakpoints with the optional usage.

from numpy import *
from pyimsl.math.constant import constant
from pyimsl.math.machine import machine
from pyimsl.math.modifiedMethodOfLines import modifiedMethodOfLines
from pyimsl.math.modifiedMethodOfLinesMgr import modifiedMethodOfLinesMgr, PDE_INITIALIZE, PDE_RESET
from pyimsl.math.writeMatrix import writeMatrix


def fcnut(npdes, x, t, u, ux, uxx, ut):
    # Define the PDE
    ut[0] = uxx[0]


def fcnbc(npdes, x, t, alpha, beta, gamp):
    alpha[0] = 0.0
    beta[0] = 1.0
    gamp[0] = 0.0


npdes = 1
nx = 10
j = 1
nstep = 10
t = [0.0]
tend = 0.0
xbreak = zeros((10), dtype='double')
y = zeros((1, 10), dtype='double')
deriv = zeros((1, 10), dtype=double)
state = []
pi = constant("pi")
arg = 9.0 * pi

# Set breakpoints and initial conditions
for i in range(0, nx):
    xbreak[i] = double(i) / double(nx - 1)
    y[0, i] = 1.0 + cos(arg * xbreak[i])
    deriv[0, i] = -arg * sin(arg * xbreak[i])

# Initialize the solver
tol = sqrt(machine(4))
modifiedMethodOfLinesMgr(PDE_INITIALIZE, state,
                         tol=tol, initialValueDerivative=deriv)

while (j <= nstep):
    tend += 0.001
    j += 1
    # Solve the problem
    modifiedMethodOfLines(t, tend, xbreak, y, state, fcnut, fcnbc)
    # Print results at every other t=tend
    if ((j % 2) != 0):
        title = "Solution at t = %5.3f" % (t[0])
        writeMatrix(title, y)
        title2 = "Derivative at t = %5.3f" % (t[0])
        writeMatrix(title2, deriv)

modifiedMethodOfLinesMgr(PDE_RESET, state)

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   -1.705e-15   -2.377e-16    4.684e-16    5.890e-15    4.064e-16
 
          7            8            9           10
  1.015e-14    1.023e-15    1.005e-14   -2.169e-19
 
                            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+00   -3.352e-15    1.974e-15    3.711e-15    3.221e-16   -3.690e-15
 
          7            8            9           10
  1.540e-15    5.639e-16    1.076e-15   -4.402e-19
 
                            Solution at t = 0.006
          1            2            3            4            5            6
      1.013        0.987        1.013        0.987        1.013        0.987
 
          7            8            9           10
      1.013        0.987        1.013        0.987
 
                           Derivative at t = 0.006
          1            2            3            4            5            6
  0.000e+00   -2.321e-15    7.869e-17    3.428e-15    8.703e-16   -3.277e-15
 
          7            8            9           10
 -3.850e-16    1.110e-15   -1.781e-13   -2.005e-12
 
                            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   -2.713e-16   -2.626e-16    3.271e-15    2.424e-16   -3.951e-15
 
          7            8            9           10
  2.005e-15    1.029e-14    1.147e-13    1.706e-12
 
                            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    5.053e-16    7.038e-16    2.135e-15   -6.885e-16   -1.424e-15
 
          7            8            9           10
 -2.229e-16   -1.293e-14   -5.296e-14   -4.409e-13

Example 4¶

In this example, consider the linear normalized hyperbolic PDE, \(u_{tt} = u_{xx}\), the “vibrating string” equation. This naturally leads to a system of first order PDEs. Define a new dependent variable \(u_t = v\). Then, \(v_t = u_{xx}\) is the second equation in the system. Take as initial data \(u(x, 0) = \sin(\pi x)\) and

\(u_t(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(\pi x) \cos(\pi t)\). Residuals are computed at the output values of t for \(0 < t \leq 2\). Output is obtained at 200 steps in increments of 0.01.

Even though the sample code modifiedMethodOfLines 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.

from __future__ import print_function
from numpy import *
from pyimsl.math.constant import constant
from pyimsl.math.machine import machine
from pyimsl.math.modifiedMethodOfLines import modifiedMethodOfLines
from pyimsl.math.modifiedMethodOfLinesMgr import modifiedMethodOfLinesMgr, PDE_INITIALIZE, PDE_RESET


def fcnut(npdes, x, t, u, ux, uxx, ut):
    # Define the PDE
    ut[0] = u[1]
    ut[1] = uxx[0]


def fcnbc(npdes, x, t, alpha, beta, 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


npdes = 2
nx = 10
j = 1
nstep = 200
t = [0.0]
tend = 0.0
xbreak = zeros((nx), dtype='double')
y = zeros((2, 10), dtype='double')
deriv = zeros((2, 10), dtype='double')
error = zeros((10), dtype=double)
state = []
pi = constant("pi")
erru = 0.0

# Set breakpoints and initial conditions
for i in range(0, nx):
    xbreak[i] = double(i) / double(nx - 1)
    y[0, i] = sin(pi * xbreak[i])
    y[1, i] = 0.0
    deriv[0, i] = pi * cos(pi * xbreak[i])
    deriv[1, i] = 0.0

# Initialize the solver
tol = sqrt(machine(4))
modifiedMethodOfLinesMgr(PDE_INITIALIZE, state,
                         tol=tol, initialValueDerivative=deriv)

while (j <= nstep):
    j += 1
    tend += 0.01
    # Solve the problem
    modifiedMethodOfLines(t, tend, xbreak, y, state, fcnut, fcnbc)
    # Look at output at steps of 0.01 and compute errors
    for i in range(0, nx):
        error[i] = y[0, i] - sin(pi * xbreak[i]) * cos(pi * tend)
        erru = max(erru, abs(error[i]))

modifiedMethodOfLinesMgr(PDE_RESET, state)
print("Maximum error in u(x,y) = %e" % (erru))

Output¶

Maximum error in u(x,y) = 1.622147e-05

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.
`IMSL_STOP_USER_FCN`	Request from user supplied function to stop algorithm. User flag = “#”.