Finite-budget coupled backward Euler for an index-one DAE

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

$$R_n(X,Z)=(X-X_{n-1}-hf(t_n,X,Z),g(t_n,X,Z))$$

and the scaled square residual

$${\mathcal{R}}_n(X,Z)=(\frac{X-X_{n-1}-hf(t_n,X,Z)}{h{\eta }_d},\frac{g(t_n,X,Z)}{{\eta }_a})\mathrm{.}$$

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.