Differential algebraic equations and consistent initialization

Let $d,m\ge 1$, let $\mathrm{\Omega }\subseteq \mathbb{R}\times {\mathbb{R}}^d\times {\mathbb{R}}^m$ be open, and let $f:\mathrm{\Omega }\to {\mathbb{R}}^d$ and $g:\mathrm{\Omega }\to {\mathbb{R}}^m$. A semi-explicit differential-algebraic equation (DAE) on $\mathrm{\Omega }$ is

The $x$-coordinates are the differential variables because their derivatives occur in the first equation. The $z$-coordinates are the algebraic variables because their derivatives do not occur. A solution on an interval $I$ is a pair of functions $x:I\to {\mathbb{R}}^d$, $z:I\to {\mathbb{R}}^m$ whose graph lies in $\mathrm{\Omega }$, with $x$ differentiable, satisfying both displayed equations at every $t\in I$. An initial pair $(x_0,z_0)$ at $t_0$ is algebraically consistent when $g(t_0,x_0,z_0)=0$. On the regular index-one branch below this is the algebraic initial condition used for reduction. For higher-index systems, algebraic consistency is only necessary: differentiated constraints may impose additional initial conditions.

At a point with $g(t_{\ast },x_{\ast },z_{\ast })=0$, the presentation is regular index one when $f,g$ are $C^1$ near the point and the square Jacobian $g_z(t_{\ast },x_{\ast },z_{\ast })=D_{z}g(t_{\ast },x_{\ast },z_{\ast })$ is invertible. The A regular index-one DAE has a local reduced ODE branch ↗ then supplies a unique local branch $z=\psi (t,x)$ through that point. Only that selected branch is used in local numerical claims. Higher-index systems and points with singular $g_z$ are outside this index-one interface.

Fix a semi-explicit presentation $x^{\mathrm{\prime }}=f(t,x,z)$, $0=g(t,x,z)$ as in Semi-explicit index-one DAEs and consistent initial data and a sufficiently differentiable consistent local solution. Differentiate the algebraic constraint along the differential equation, substituting $x^{\mathrm{\prime }}=f(t,x,z)$ and any previously obtained derivative relations at each stage. The differentiation index in this convention is the least positive integer $k$ for which the original constraint and its first $k$ differentiated equations locally determine $z$ and $z^{\mathrm{\prime }}$, with the equation for the highest derivative of $z$ having an invertible square Jacobian on a neighborhood of the solution. If no such finite $k$ exists, this local index is undefined. This convention counts through the first derivative of the algebraic variables and is local to the fixed presentation; it is not a universal index theory for arbitrary implicit systems.

If $g_z$ is invertible, one differentiation gives

which determines $z^{\mathrm{\prime }}$, so the index is one. If $g_z$ is singular, the index must be checked from the actual differentiated equations; no index-one conclusion follows from this definition alone.

Let $d,m\ge 1$, let $\mathrm{\Omega }\subseteq \mathbb{R}\times {\mathbb{R}}^d\times {\mathbb{R}}^m$ be open, and let $f:\mathrm{\Omega }\to {\mathbb{R}}^d$, $g:\mathrm{\Omega }\to {\mathbb{R}}^m$ be $C^k$ for $k\ge 1$. Suppose $g(t_{\ast },x_{\ast },z_{\ast })=0$ and $g_z(t_{\ast },x_{\ast },z_{\ast })$ is invertible. Then there are neighborhoods $U$ of $(t_{\ast },x_{\ast })$ and $V$ of $z_{\ast }$, and a unique $C^k$ map $\psi :U\to V$, such that, for $(t,x)\in U$ and $z\in V$,

Its derivatives satisfy

with the derivatives of $g$ evaluated at $(t,x,\psi (t,x))$. The reduced vector field $F(t,x)=f(t,x,\psi (t,x))$ is $C^k$ on $U$. Thus solutions of the DAE that stay in this selected chart are exactly solutions of $x^{\mathrm{\prime }}=F(t,x)$ paired with $z=\psi (t,x)$.

Let $K$ be a compact set of base points $(t,x)$ in a selected branch chart from A regular index-one DAE has a local reduced ODE branch, and put $z_{\ast }=\psi (t,x)$. Suppose $g(t,x,\cdot )$ is $C^2$ on $B(z_{\ast },\rho )$ for every $(t,x)\in K$. Let $G_{\ast }(t,x)=g_z(t,x,\psi (t,x))$. Assume common constants $\gamma >0$, $\mathrm{\Lambda }\ge 0$, and $\rho >0$ satisfy

for every $(t,x)\in K$ and $z\in B(\psi (t,x),\rho )$; if $\mathrm{\Lambda }>0$, assume $\rho \le \gamma \mathrm{/}(\sqrt{m}\mathrm{\Lambda })$. Then, uniformly for all such base points and candidates,

The estimate is local: it applies only to candidates in the specified branch ball and does not turn a residual into a global branch or uniqueness certificate.

Let $f,g$ define a regular index-one DAE and let $U\times V$ be a selected local branch chart from A regular index-one DAE has a local reduced ODE branch, so $g(t,x,z)=0$ is equivalent to $z=\psi (t,x)$ there. Input $t_0,x_0$, a seed $z_{\mathrm{s}\mathrm{e}\mathrm{e}\mathrm{d}}\in V$, an algebraic residual tolerance ${\epsilon }_a\ge 0$, and a finite Newton-correction budget $N\in {\mathbb{N}}_0$. If $(t_0,x_0)\mathrm{\notin }U$, $g(t_0,x_0,\cdot )$ is undefined at the seed, or dimensions do not match, return invalid-input.

Call Finite-budget Newton for a square vector residual on the square residual $G(z)=g(t_0,x_0,z)$, with initial value $z_{\mathrm{s}\mathrm{e}\mathrm{e}\mathrm{d}}$, tolerance ${\epsilon }_a$, and budget $N$, restricting evaluations to the chart's open $z$-domain. Write $z_0$ for its returned last valid iterate. If it reports root, reevaluate $G(z_0)$ once: if that residual is exactly zero, return $(x_0,z_0,\mathtt{\text{exact-consistent}})$; otherwise return $(x_0,z_0,\mathtt{\text{approximate-consistent}})$ and its norm. Map singular-jacobian, invalid-trial, and iteration-budget-exhausted to the corresponding initialization statuses and return the last valid iterate.

The procedure makes no arbitrary-start convergence claim. approximate-consistent records a residual bound; it is not relabeled as exact consistency. Under the index-one branch hypotheses an exact consistent value exists locally, by A regular index-one DAE has a local reduced ODE branch, but the finite Newton budget may fail to find it.

The solve has $m$ algebraic unknowns. Record residual queries, Jacobian assemblies and correction solves separately. Their costs depend on the chosen backend.

Input $f,g$, a constant grid $t_n=t_0+nh$ with $h>0$ and $N_t\in \mathbb{N}$, an initial pair $(X_0,Z_0)$, and a selected open tube $\mathcal{T}\subset \mathrm{\Omega }$ inside the local branch chart of A regular index-one DAE has a local reduced ODE branch, where $g(t,x,z)=0$ is equivalent to $z=\psi (t,x)$. For each time $t$, write ${\mathcal{T}}_t=\{(X,Z):(t,X,Z)\in \mathcal{T}\}$ for the open state slice used as the Newton evaluation domain. Also input positive residual tolerances ${\eta }_d,{\eta }_a$, a finite Newton-correction budget $B\in {\mathbb{N}}_0$, and an initial algebraic-residual bound ${\eta }_{a,0}\ge 0$. Return invalid-grid if $h\le 0$ or $N_t=0$, and invalid-input for incompatible dimensions, data outside the function domains, or $(X_0,Z_0)\mathrm{\notin }{\mathcal{T}}_{t_0}$. Compute the starting algebraic residual once; return inconsistent-initial-data if $\mathrm{\parallel }g(t_0,X_0,Z_0){\mathrm{\parallel }}_2>{\eta }_{a,0}$.

For $n=1,\dots ,N_t$, form the coupled equation residual

and the scaled square residual

Use $(X_{n-1},Z_{n-1})$ as the Newton seed. Before calling Newton, if $(X_{n-1},Z_{n-1})\mathrm{\notin }{\mathcal{T}}_{t_n}$, return invalid-trial with the already computed prefix. Otherwise call Finite-budget Newton for a square vector residual on ${\mathcal{R}}_n$ with residual tolerance $1$, budget $B$, and ${\mathcal{T}}_{t_n}$ as its evaluation domain. If the returned status is root and the candidate lies in ${\mathcal{T}}_{t_n}$, reevaluate $R_n$ once to retain both unscaled residuals, then append the candidate. Otherwise return the corresponding singular-Jacobian, invalid-trial, or iteration-budget-exhausted status and the already computed prefix. The input functions are defined on the selected tube, so no separate evaluation-failure branch is needed. If every step is accepted, return completed and the state sequence.

For recording, define the derivative residual $r_n^d=(X_n-X_{n-1})\mathrm{/}h-f(t_n,X_n,Z_n)$ and the algebraic residual $r_n^a=g(t_n,X_n,Z_n)$. root means only that the scaled residual tolerance was met; it does not assert an exact zero, uniqueness, Newton convergence from an arbitrary seed, or a discretization error bound.

Use the coupled $(d+m)$-dimensional residual and its current block Jacobian. Record residual queries, Jacobian assemblies and correction solves separately; backend costs are not specified.

Input $f,g$, a constant grid $t_n=t_0+nh$ with $h>0$ and $N_t\ge 2$, two starting pairs $(X_0,Z_0),(X_1,Z_1)$, and a selected open tube $\mathcal{T}\subset \mathrm{\Omega }$ inside the local branch chart of A regular index-one DAE has a local reduced ODE branch, where $g(t,x,z)=0$ is equivalent to $z=\psi (t,x)$. For each time $t$, write ${\mathcal{T}}_t=\{(X,Z):(t,X,Z)\in \mathcal{T}\}$ for the open state slice used as the Newton evaluation domain. Also input positive residual tolerances ${\eta }_d,{\eta }_a$, a finite Newton-correction budget $B\in {\mathbb{N}}_0$, and two starting algebraic-residual bounds ${\eta }_{a,0},{\eta }_{a,1}\ge 0$. Return invalid-grid for a nonpositive step or fewer than two starting states, and invalid-input for incompatible dimensions, data outside the function domains, or $(X_j,Z_j)\mathrm{\notin }{\mathcal{T}}_{t_j}$ for either $j=0,1$. Compute $r_j^a=g(t_j,X_j,Z_j)$ for both $j=0,1$; if $\mathrm{\parallel }{r}_j^a{\mathrm{\parallel }}_2>{\eta }_{a,j}$ for either start, return inconsistent-initial-data.

For $n=2,\dots ,N_t$, form

and

Use $(2X_{n-1}-X_{n-2},Z_{n-1})$ as the Newton seed. Before calling Newton, if $(2X_{n-1}-X_{n-2},Z_{n-1})\mathrm{\notin }{\mathcal{T}}_{t_n}$, return invalid-trial with the computed prefix. Otherwise call Finite-budget Newton for a square vector residual on ${\mathcal{R}}_n$ with tolerance $1$, budget $B$, and ${\mathcal{T}}_{t_n}$ as the evaluation domain. Append a candidate only when Newton reports root and the candidate lies in ${\mathcal{T}}_{t_n}$; reevaluate $R_n$ once to record its unscaled residual blocks. Otherwise return the corresponding singular-Jacobian, invalid-trial, or iteration-budget-exhausted status with the computed prefix. The input functions are defined on the selected tube, so no separate evaluation-failure branch is needed. If all steps are accepted, return completed and the state sequence $((X_0,Z_0),\dots ,(X_{N_t},Z_{N_t}))$.

For an accepted step, the constant-step BDF2 derivative residual is

and the algebraic residual is $r_n^a=g(t_n,X_n,Z_n)$. The procedure makes no claim that Newton finds a root or that an accepted state is the exact selected branch root.

Use the coupled $(d+m)$-dimensional residual and its current block Jacobian. Record residual queries, Jacobian assemblies and correction solves separately; backend costs are not specified.

Let a compact exact solution segment lie in a common index-one branch chart with $f,g\in C^1$ and $g_z$ invertible on the graph. For each point $(t,x,\psi (t,x))$ on that segment, the coupled backward-Euler root equations

have a unique root $(X,Z)$ in a selected neighborhood of $(x,\psi (t,x))$ whenever the step data $(h,t,\xi )$ are sufficiently close to $(0,t,x)$. On a compact subsegment contained in one chart with positive margin, the neighborhoods and a positive step bound can be chosen uniformly after restricting the input states to that neighborhood.

Let a compact exact solution segment lie in a common index-one branch chart with $f,g\in C^1$ and $g_z$ invertible on the graph. For each point $(t,x,\psi (t,x))$, the coupled constant-step BDF2 root equations

have a unique root $(X,Z)$ in a selected neighborhood of $(x,\psi (t,x))$ whenever $(h,t,\xi ,\zeta )$ is sufficiently close to $(0,t,x,x)$. On a compact subsegment in one chart with positive margin, one may choose a common positive step bound after restricting the two input states to neighborhoods of the reference graph.

In the exact-arithmetic model with returning residual, Jacobian, domain-test and correction primitives, suppose the initialization and stepper algorithms in Finite-budget index-one DAE consistent initialization, Finite-budget coupled backward Euler for an index-one DAE, and Finite-budget coupled constant-step BDF2 for an index-one DAE are run with finite grid length and finite Newton budget $B$. Count explicit residual queries separately from work internal to Jacobian assembly and linear-system backends. Each call returns a status after finitely many such queries. A Newton call makes at most $B$ correction calls and $B+1$ explicit residual evaluations. Initialization makes at most one additional algebraic-residual evaluation after a root status to record the returned residual norm, so it uses at most $B+2$ evaluations. A backward-Euler call with $S$ attempted steps uses at most $1+S(B+2)$ residual evaluations, including its one starting algebraic-residual check; a BDF2 call with $S$ attempted steps uses at most $2+S(B+2)$, including its two starting algebraic-residual checks. These bounds include calls that return early with a failure status. Every successful correction solves one current Jacobian system, of dimension $m$ in initialization and $d+m$ in either stepper; assembly and solve costs depend on the backend.

For either stepper, every appended candidate satisfies $\mathrm{\parallel }{r}_n^d{\mathrm{\parallel }}_2\le {\eta }_d$ and $\mathrm{\parallel }{r}_n^a{\mathrm{\parallel }}_2\le {\eta }_a$ for its defined derivative and algebraic residuals. These are conditional accepted-residual bounds. No residual tolerance is guaranteed to be attained, and these bounds assert neither an exact zero, uniqueness nor a discretization error estimate.

For the semi-explicit DAE with $x\in {\mathbb{R}}^d$, $z\in {\mathbb{R}}^m$, and $d,m\ge 1$, let $x^{\mathrm{\prime }}=F(t,x)$ be the reduced ODE from A regular index-one DAE has a local reduced ODE branch on $[0,T]$, with exact solution $\bar{x}\in C^2$. Suppose all exact and numerical states lie in a common compact branch tube, $F$ is Lipschitz in $x$ there with constant $L$, $\psi$ is Lipschitz in $x$ with constant $L_{\psi }$, and $f$ is Lipschitz in $z$ with constant $L_z$. Assume the residual lemma's hypotheses hold uniformly on this tube: there are common $\gamma >0$, $\mathrm{\Lambda }\ge 0$, and $\rho >0$ for its lower singular-value and Jacobian-Lipschitz bounds at every base point, with $\rho \le \gamma \mathrm{/}(\sqrt{m}\mathrm{\Lambda })$ when $\mathrm{\Lambda }>0$; all numerical algebraic candidates lie in the corresponding balls $B(\psi (t_n,X_n),\rho )$. Use a uniform grid $t_n=nh$, $h=T\mathrm{/}N$, $T>0$, and $N\ge 1$. Assume the starting pair also lies in its branch ball and

For each $1\le n\le N$, suppose $(X_n,Z_n)$ lies in the branch ball and its coupled backward-Euler residuals satisfy

Then a constant $C$, independent of $h$, satisfies

for all sufficiently small $h$. This is a conditional accuracy result for accepted residuals; the finite-budget solver may return a failure status.

For the semi-explicit DAE with $x\in {\mathbb{R}}^d$, $z\in {\mathbb{R}}^m$, and $d,m\ge 1$, let $x^{\mathrm{\prime }}=F(t,x)$ be the reduced ODE from A regular index-one DAE has a local reduced ODE branch on $[0,T]$, with exact solution $\bar{x}\in C^3$. Assume $F$ is Lipschitz in $x$ on a common compact tube with constant $L$, $\psi$ is Lipschitz in $x$ there with constant $L_{\psi }$, $f$ is Lipschitz in $z$ with constant $L_z$, and the residual lemma's hypotheses hold uniformly on this tube: common $\gamma >0$, $\mathrm{\Lambda }\ge 0$, and $\rho >0$ in its lower singular-value and Jacobian-Lipschitz bounds, with $\rho \le \gamma \mathrm{/}(\sqrt{m}\mathrm{\Lambda })$ when $\mathrm{\Lambda }>0$. Every numerical algebraic candidate lies in its corresponding branch ball. Use a constant grid $t_n=nh$, $h=T\mathrm{/}N$, $T>0$, and $N\ge 2$. Suppose the two starting pairs satisfy

For each $2\le n\le N$, suppose the constant-step BDF2 derivative and algebraic residuals satisfy

Then, for all sufficiently small $h$, a constant $C$ independent of $h$ satisfies

The starter hypothesis holds, in particular, for an exact consistent $X_0,Z_0$ and one backward-Euler startup from $X_0$ with derivative residual $O(h)$ and algebraic residual $O(h^2)$. The error theorem is conditional on the residual bounds and selected branch; it does not assert that a finite Newton budget meets them.