Finite-budget coupled constant-step BDF2 for an index-one DAE

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

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

and

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

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

$$r_n^d=\frac{3X_n-4X_{n-1}+X_{n-2}}{2h}-f(t_n,X_n,Z_n),$$

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.