Introduction to pde1dMg

The section describes an algorithm and a corresponding integrator function pde1dMg for solving a system of partial differential equations

Equation 1

$u_t \equiv \frac{\partial u}{\partial t} = f (u, x, t), x_L < x < x_R, t > t_0$

This software is a one-dimensional differential equation solver. It requires the user to provide initial and boundary conditions in addition to a function for the evaluation of $u_t$. The integration method is noteworthy due to the maintenance of grid lines in the space variable, $x$. Details for choosing new grid lines are given in Blom and Zegeling, (1994). The class of problems solved with pde1dMg is expressed by Equation 1 and given in more detail by:

Equation 2

$$\begin{split}\begin{array}{c} \displaystyle\sum_{k=1}^{\textit{NPDE}} C_{j,k}\left(x, t, u, u_x\right)\frac{{\partial}u^k}{{\partial}t} = x^{-m} \frac{\partial}{{\partial}x}\left(x^mR_j\left(x, t, u, u_x\right)\right) - Q_j\left(x, t, u, u_x\right), \\ j = 1, \ldots, \textit{NPDE}, \phantom{...} x_L < x < x_R, \phantom{...} t > t_0, \phantom{...} m \in { 0, 1, 2 } \end{array}\end{split}$$

The vector $u \equiv \left[ u^1, \ldots, u^{\mathit{NPDE}} \right]^T$ is the solution. The integer value $\mathit{NPDE} \ge 1$ is the number of differential equations. The functions $R_j$ and $Q_j$ can be regarded, in special cases, as flux and source terms. The functions $u, C_{j,k}, R_j, Q_j$ are expected to be continuous. Allowed values for the integer $m$ are any of $m = 0, 1, 2$. These are respectively for problems in Cartesian, cylindrical or polar, and spherical coordinates. In the two cases with $m > 0$, the interval $\left[ x_L, x_R \right]$ must not contain $x = 0$ as an interior point.

The boundary conditions have the master equation form

Equation 3

$$\begin{split}\begin{array}{l} \beta_j(x,t) R_j\left(x,t,u,u_x\right) = \gamma_j\left(x,t,u,u_x\right), \\ \text{at } x = x_L \text{ and } x = x_R, j = 1, \ldots \mathit{NPDE} \end{array}\end{split}$$

In the boundary conditions the functions $\beta_j$ and $\gamma_j$ are continuous. In the two cases with $m > 0$, with an endpoint of $\left[ x_L, x_R \right]$ at 0, the finite value of the solution at $x = 0$ must be ensured. This requires the specification of the solution at $x = 0$, or it implies that $R_j | _{x = x_L} = 0$ or $R_j | _{x = x_R} = 0$. The initial values satisfy $u \left( x, t_0 \right) = u_0 (x), \phantom{...} x \in \left[ x_L, x_R \right]$ , where $u_0$ is a piece-wise continuous vector function of $x$ with $\mathit{NPDE}$ components.

The user must pose the problem so that mathematical definitions are known for the functions

$$C_{j,k}, R_j, Q_j, \beta_j, \gamma_j \text{ and } u_0$$

These functions are provided to the function pde1dMg in the form of two user-supplied functions. This form of the usage interface is explained below and illustrated with several examples. $u_0$ can be supplied as the input argument u or by an optional user-supplied function. Users comfortable with the description of this algorithm may skip directly to the Examples section.

Description Summary

Equation 1 is approximated at $N = \mathit{ngrids}$ time-dependent grid values $x_L = x_0 < x_1 < \ldots < x_i(t) < \ldots < x_{N+1} = x_R$. Using the total differential $\tfrac{du}{dt} = u_t + u_x \tfrac{dx}{dt}$ transforms the differential equation to the form

$$u_t = \frac{du}{dt} - u_x\frac{dx}{dt} = f(u,x,t), \phantom{...} x_L < x < x_R$$

Using central divided differences for the factor $u_x$ leads to the system of ordinary differential equations in implicit form

$$\frac{dU_{\mathrm{i}}}{dt} - \frac {\left(U_{\mathrm{i}+1} - U_{\mathrm{i}-1}\right)} {\left(x_{\mathrm{i}+1} - x_{\mathrm{i}-1}\right)} \frac{dx_{\mathrm{i}}}{dt} = F_{\mathrm{i}}, \phantom{...} t > t_0, \phantom{...} i = 1, \cdots, N$$

The terms $U_i, F_i$ respectively represent the approximate solution to the partial differential equation and the value of $f (u, x, t)$ at the point $(u, x, t) = \left( u_i, x_i(t), t \right)$. The truncation error from this approximation is second-order in the space variable $x$. The above ordinary differential equations are underdetermined, so additional equations are added for determining the time-dependent grid points. These additional equations contain parameters that can be adjusted by the user. Often it will be necessary to modify these parameters to solve a difficult problem. For this purpose the following quantities are needed:

$$\begin{split}\begin{array}{l} \mathit{\Delta}x_i = x_{i+1} - x_i, \phantom{...} n_i = \mathit{\Delta}x_i^{-1} \\ \mu_i = n_i - \kappa(\kappa+1)\left(n_{i+1}-2n_i + n_{i-1}\right), \phantom{...} 0 \leq i \leq N \\ n_{-1} \equiv n_0, \phantom{...} n_{N+1} \equiv n_N \end{array}\end{split}$$

The values $n_i$ are the so-called point concentration of the grid. The parameter $\kappa \geq 0$ denotes a spatial smoothing value. Now the grid points are defined implicitly so that

$$\frac{\mu_{i-1} + \tau \tfrac{du_{i-1}}{dt}}{M_{i-1}} = \frac{\mu_i + \tau \tfrac{du_i}{dt}}{M_i}, \phantom{...} 1 \leq i \leq N$$

The parameter $\tau \geq 0$ denotes a time-smoothing value. If the value $\tau$ is chosen to be large, this results in a fixed spatial grid. Increasing $\tau$ from its default value avoids the error condition where grid lines cross. The divisors are defined by

$$M_i^2 = \alpha + \mathit{NPDE}^{-1} \sum_{j=1}^{\mathit{NPDE}} \frac{\left(U_{i+1}^j - U_i^j\right)^2}{\left(\mathit{\Delta}x_i\right)^2}$$

The value $\kappa$ determines the level of clustering or spatial smoothing of the grid points. Decreasing $\kappa$ from its default values also decreases the amount of spatial smoothing. The parameters $M_i$ approximate arc length and help determine the shape of the grid or $x_i$ distribution. The parameter $\tau$ prevents the grid movement from adjusting immediately to new values of the $M_i$ , thereby avoiding oscillations in the grid that cause large relative errors in the solution. This is important when applied to solutions with steep gradients.

The discrete form of the differential equation and the smoothing equations are combined to yield the implicit system of differential equations

$$\begin{split}\begin{array}{l} A(Y)\tfrac{dY}{dt} = L(Y), \\ Y = \left[ U_1^1, \ldots, U_1^{\mathit{NPDE}}, x_1, \ldots \right]^T \end{array}\end{split}$$

This is usually a stiff differential-algebraic system. It is solved using the integrator differentialAlgebraicEqs. If differentialAlgebraicEqs is needed during the evaluations of the differential equations or boundary conditions, it must be done in a separate thread to avoid possible problems with pde1dMg’s internal use of differentialAlgebraicEqs. The only options for differentialAlgebraicEqs set by pde1dMg are the Maximum BDF Order, and the absolute and relative error values, documented as maxBdfOrder, and atolRtolScalars.