6 Construction and Analysis of RKM for nonlinear DAEs

Now we consider RKM for nonlinear DAEs. We start with a DAE in semi explicit strangeness-free form and give general results on how to write down a general RKM for it and how to analyse the global error. Then, we consider general strangeness-free nonlinear DAEs and show that a certain class of RKM applies well – namely those that can be constructed by collocation with Lagrange polynomials over the Radau, Lobatto, or Gauss quadrature points.

6.1 General RKM for Semi-Explicit Strangeness-free DAEs

A semi explicit strangeness-free DAE is of the form \[\begin{align} \dot x &= f(t, x, y) \tag{6.1} \\ 0 &= g(t, x, y) \tag{6.2} \end{align}\] with the Jacobian of \(g\) with respect to \(y\), i.e. \[ \partial_y\otimes g(t, x(t), y(t)) =: g_y(t, x(t), y(t)), \] being invertible for all \(t\) along the solution \((x,y)\).

Some observations:

  • this system is strangeness-free
  • under certain assumptions, any DAE can be brought into this form
  • in the linear case \(E\dot z = Az +f\), with \(z=(x,y)\), the assumptions basically mean that \[ E = \begin{bmatrix} I & 0 \\ 0 & 0 \end{bmatrix} \quad\text{and}\quad A = \begin{bmatrix} * & * \\ * & A_{22} \end{bmatrix}, \] with \(A_{22}(t)\) invertible for all \(t\).
  • The condition \(g_y\) invertible means that, locally, one could consider \[ \dot x = f(t, x, R(t,x)), \quad\text{with $R$ such that}\quad y=R(t,x). \] However, this is not practical for numerical purposes.

The general strategy to get a suitable formulation of a time discretization of system (6.1)-(6.2) by any RKM is to consider the perturbed version \[\begin{align*} \dot x = f(t, x, y), \\ \varepsilon \dot y = g(t, x, y), \end{align*}\] which is an ODE, formulate the RKM, and then let \(\varepsilon \to 0\). In the Hairer/Wanner Book, this approach is called *\(\varepsilon\)-embedding.

This is, consider

Table 6.1: RKM applied to semi-explicit DAEs
\(x_{i+1} = x_i + h\sum_{j=1}^s\beta_j \dot X_{ij}\), \(y_{i+1} = y_i + h\sum_{j=1}^s\beta_j \dot Y_{ij}\),
\(\dot X_{ij} = f(t_i+\gamma_jh, X_{ij}, Y_{ij})\), \(\varepsilon \dot Y_{ij} = g(t_i+\gamma_j h, X_{ij}, Y_{ij})\), \(j=1,2,\cdots,s, \quad \quad (*)\)
\(X_{ij} = x_i + h\sum_{\ell=1}^s\alpha_{j\ell}\dot X_{i\ell}\), \(\phantom{\varepsilon}Y_{ij} = y_i + h\sum_{\ell=1}^s\alpha_{j\ell}\dot Y_{i\ell}\), \(j=1,2,\cdots,s,\)

i.e., the RKM applied to an ODE in the variables \((x,y)\), and replace \((*)\) by

\[ \dot X_{ij} = f(t_i+\gamma_jh, X_{ij}, Y_{ij}),\quad 0 = g(t_i+\gamma_j h, X_{ij}, Y_{ij}), \quad j=1,2,\cdots,s. \]

Theorem 6.1 (Kunkel/Mehrmann Thm. 5.16) Consider a semi-explicit, strangeness-free DAE as in (6.1)-(6.2) with a consistent initial value \((x_0, y_0)\). The time-discretization by a RKM,

  • with \(\mathcal A\) invertible and \(\rho:=1-\beta^T\mathcal A^{-1}e\),
  • applied as in Table 6.1 with \(\varepsilon=0\),
  • that is convergent of order \(p\) for ODEs
  • and fulfills the Butcher condition \(C(q)\) with \(q\geq p+1\)

leads to an global error that behaves like \[ \|\mathfrak X(t_N) - \mathfrak X_N\| = \mathcal O(h^k), \]


  • \(k=p\), if \(\rho=0\),
  • \(k=\min\{p, q+1\}\), if \(-1\leq \rho < 1\)
  • \(k=\min\{p, q-1\}\), if \(\rho =1\).

If \(|\rho|>1\), then the RKM – applied to (6.1)(6.2) – does not converge.

Some words on the conditions on \(p\), \(q\), and \(\rho\):

  • For stiffly accurate methods, \(\beta^T \mathcal A^{-1}e=1\) and, thus, \(\rho=0\) \(\rightarrow\) no order reduction for strangeness free or index-1 systems
  • For the implicit midpoint rule also known as the 1-stage Gauss method: \[ \begin{array}{c|c} \frac 12 & \frac 12 \\ \hline & 1 \end{array} \]

    • the convergence order for ODEs is \(p=2\)
    • but \(1-\beta^T \mathcal A^{-1}e = 1- 1\cdot {\bigl(\frac 12\bigr)}^{-1} 1 = -1\), so that \(\min\{p-1, q\} = k \leq 1\), depending on \(q\).
    • in fact \(k=1\) as \[ C(q): \quad \sum_{\ell=1}^s\alpha_{j\ell}\gamma_\ell^{\bar k-1}=\frac{1}{\bar k} \gamma_j^{\bar k}, \quad {\bar k}=1,\cdots,q, \quad j=1,\cdots,s \] in the present case of \(s=1\), \(\alpha_{11}=\gamma_1=\frac 12\) is fulfilled for \({\bar k}=1: \quad \frac 12 = \frac 12\)
    • it is not relevant here, but for \({\bar k}=2:\quad \frac 12 \cdot \frac 12 \neq \frac 12 \cdot \frac 14\)

6.2 Collocation RKM for Implicit Strangeness-free DAEs

The general form of a strangeness-free DAE is given as \[\begin{align} \hat F_1(t,x,\dot x) &= 0 \tag{6.3}\\ \hat F_2(t,x) &= 0 \tag{6.4} \end{align}\] where the strangeness-free or index-1 assumption is encoded in the existence of implicit functions \(\mathcal L\), \(\mathcal R\) such that, with \(x=(x_1,x_2)\), the implicit DAE (6.3)(6.4) is equivalent to the semi-explicit DAE \[\begin{align*} \dot x_1 &= \mathcal L(t,x_1,x_2) \\ 0 &= \mathcal R(t,x_1) -x_2 \end{align*}\]

In what follows we show that a collocation approach coincides with certain RKM discretizations so that the convergence analysis of the RKM can be done via approximation theory.

Regression (Collocation): – If one looks for a function \(x\colon [0,1] \to \mathbb R^{}\) that fulfills \(F(x(t))=0\) for all \(t\in[0,1]\), one may interpolate \(x\) by, say, a polynomial \(x_p(t) = \sum_{\ell=0}^kx_\ell t^\ell\) and determine the \(k+1\) coefficients \(x_\ell\) via the solution of the system of (nonlinear) equations \(F(x_p(t_\ell))=0\), \(\ell=0,1,\dotsc,k\), where the \(t_\ell\in[0,1]\) are the \(k+1\) collocation points.

Concretely, we parametrize \(s\) collocation points via \[\begin{equation} 0 = \gamma_0 < \gamma_1 <\gamma_2< \dotsc < \gamma_s=1 \tag{6.5} \end{equation}\] and define two sets of Lagrange polynomials \[ L_\ell(\xi) = \prod_{j=0,j\neq \ell}^s \frac{\xi-\gamma_j}{\gamma_\ell-\gamma_j} \quad\text{and}\quad \tilde L_\ell(\xi) = \prod_{m=1,m\neq \ell}^s \frac{\xi-\gamma_m}{\gamma_\ell-\gamma_m}, \] with \(\ell\in\{0,1,\dotsc,s\}\).

Let \(\mathbb P_k\) be the space of polynomials of degree \(\leq k-1\). We define the collocation polynomial \(x_\pi \in \mathbb P_{s+1}\) via \[\begin{equation} x_\pi (t) = \sum_{\ell=0}^s X_{i\ell}L_\ell\bigl(\frac{t-t_i}{h}\bigr) \tag{6.6} \end{equation}\] designed to compute the stage values \(X_{i\ell}\), where \(X_{i0}=x_i\) is already given.

The stage derivatives are then defined as \[\begin{equation} \dot X_{ij} = \dot x_\pi(t_i+\gamma_jh) = \frac 1h \sum_{\ell=0}^sX_{i\ell}\dot L_\ell(\gamma_j). \tag{6.7} \end{equation}\]

To obtain \(x_{i+1}=x_\pi(t_{i+1})=X_{is}\), we require the polynomial to satisfy the DAE (6.3)(6.4) at the collocation points \(t_{ij}=t_i+\gamma_jh\), that is

\[\begin{equation} \hat F_1(t_i+\gamma_jh,X_{ij},\dot X_{ij}) = 0, \quad \hat F_2(t_i+\gamma_jh,X_{ij}) = 0, \quad j=1,\dotsc,s. \phantom{F_1} \tag{6.8} \end{equation}\]

Now we show that this collocation defines a RKM discretization of (6.3)(6.4).

Since \(\tilde L_\ell \in \mathbb P_s\) for \(\ell=1,\ldots,s\) it holds that \[ P_\ell(\sigma):=\int_0^\sigma \tilde L_\ell (\xi)d\xi \in \mathbb P_{s+1} \] that is, by Lagrange interpolation, it can be written as \[ P_\ell(\sigma) = \sum_{j=0}^s P_\ell(\gamma_j)L_j(\sigma). \]

If we differentiate \(P_\ell\), we get \[ \dot P_\ell(\sigma) = \sum_{j=0}^s P_\ell(\gamma_j)\dot L_j(\sigma) = \sum_{j=0}^s \int_0^{\gamma_j} \tilde L_\ell (\xi)d\xi \dot L_j(\sigma)=: \sum_{j=0}^s \alpha_{j\ell} \dot L_j(\sigma) \] where define simply define \[ \alpha_{j\ell} = \int_0^{\gamma_j} \tilde L_\ell (\xi)d\xi. \] Note that \(\alpha_{0\ell}=0\) such that summation in \(\dot P_\ell(\sigma)\) starts at \(j=1\) instead of \(j=0\). Moreover, by definition of \(P_\ell\) (and the fundamental theorem of calculus), it holds that \[ \dot P_\ell(\sigma) = \tilde L_\ell(\sigma), \] which gives that \(\dot P_\ell(\gamma_m) = \delta_{\ell m}\) that is \[ \dot P_\ell(\gamma_m) = \sum_{j=1}^s\alpha_{j\ell}\dot L_j(\gamma_m) = \begin{cases} 1, &\quad \text{if }\ell =m \\ 0, &\quad \text{otherwise} \end{cases} \] for \(\ell, m=1,\dotsc,s\).

Accordingly, if we define \(\mathcal A := \bigl[\alpha_{j\ell}\bigr]_{j,\ell=1,\dotsc,s} \in \mathbb R^{s,s}\) and \[ V:=\bigl[v_{mj}\bigr]_{m,j=1,\dotsc,s} = \bigl[ \dot L_j(\gamma_m) \bigr]_{m,j=1,\dotsc,s} \in \mathbb R^{s,s} , \] it follows that \(V=\mathcal A^{-1}\).

Moreover, since, \[ \sum_{j=0}^s L_j(\sigma) \equiv 1, \quad\text{so that }\quad\sum_{j=0}^s \dot L_j(\sigma) \equiv 0, \] we have that \[ \sum_{j=0}^s \dot L_j(\gamma_m) =0= \sum_{j=0}^s v_{mj} \] and, thus, \[ v_{m0} = -\sum_{j=1}^s \dot L_j(\gamma_m) = -e_m^TVe. \] With these relations we rewrite (6.7) as \[ h\dot X_{im} = \sum_{\ell=0}^sX_{i\ell}\dot L_\ell(\gamma_m) = v_{m0}x_i + \sum_{\ell=1}^sv_{m\ell}X_{i\ell}. \] and \(h\sum_{m=1}^s\alpha_{\ell m} \dot X_{im}\) as \[\begin{align} h\sum_{m=1}^s \alpha_{\ell m}\dot X_{im} &= \sum_{m=1}^s \alpha_{\ell m}v_{m0}x_i + \sum_{j,m=1}^s \alpha_{\ell m}v_{mj}X_{ij} \notag \\ &= -e_\ell^T \mathcal AV e x_i + \sum_{j=1}^se_\ell^T \mathcal AVe_jX_{ij} \tag{6.9}\\ &= -x_i + X_{i\ell}, \notag \end{align}\] which, together with (6.8), indeed defines a RKM.

Some remarks:

  • the preceding derivation shows that the collocation (6.6) and (6.8) is equivalent to the RKM scheme (6.9) and (6.8)
  • convergence of these schemes applied to (6.3)(6.4) is proven in Kunkel/Mehrmann Theorem 5.17
  • with fixing \(\gamma_s=1\), the obtained RKM is stiffly accurate
  • the remaining \(s-1\) \(\gamma\)s can be chosen to get optimal convergence rates \(\rightarrow\) RadauIIa schemes
  • if also \(\gamma_s\) is chosen optimal in terms of convergence, the Gauss schemes are obtained