Compatible matrix evolution and time stepping

In positive dimension, if $M_{e}x=0$, then $x^T{M}_{e}x=0$, so positive definiteness gives $x=0$; the same argument applies to $M_b$. In dimension zero the unique square matrix is the identity map on the singleton space and is invertible. Thus both matrices are injective and hence invertible in their finite-dimensional domains, so their inverse operations used by consumers are defined in every allowed dimension.

For $u\in {\mathbb{R}}^{n_b}$ and $v\in \ker C$, $⟨C^{T}u,v⟩=u^{T}Cv=0$, hence $\mathrm{im}{C}^T\subseteq (\ker C{)}^{\perp }$. Also $z\in (\mathrm{im}{C}^T{)}^{\perp }$ exactly when $u^{T}Cz=0$ for every $u\in {\mathbb{R}}^{n_b}$, which is exactly $Cz=0$; therefore $(\mathrm{im}{C}^T{)}^{\perp }=\ker C$. By [F1], decompose any $w\in (\ker C{)}^{\perp }$ as $w=u+v$ with $u\in \mathrm{im}{C}^T$ and $v\in (\mathrm{im}{C}^T{)}^{\perp }=\ker C$. Both $w$ and $u$ are orthogonal to $v$, so $\parallel v{\parallel }^2=⟨w-u,v⟩=0$ and $v=0$. Hence $(\ker C{)}^{\perp }=\mathrm{im}{C}^T$.

Assume $\mathrm{im}G=\ker C$. If $q=0$, then $G^T{M}_{e}e=0$, so $M_{e}e$ is orthogonal to every column of $G$ and hence to $\ker C$. By step 1.2, $M_{e}e\in \mathrm{im}{C}^T$.

Conversely, if $M_{e}e\in \mathrm{im}{C}^T$, then step 1.2 gives $M_{e}e\in (\ker C{)}^{\perp }$. Every column of $G$ lies in $\ker C$, so $G^T{M}_{e}e=0$ and therefore $q=0$. The two implications prove the stated characterization.

Put $z=(e,b)\in {\mathbb{R}}^{n_e+n_b}$ and define $A:=\left(\begin{array}{cc}-M_e^{-1}S& M_e^{-1}{C}^T{M}_b\\ -C& 0\end{array}\right)$ and $f(t):=\left(\begin{array}{c}-M_e^{-1}j(t)\\ 0\end{array}\right)$. The inverses exist by positive definiteness in [F1]: if $M_{e}x=0$, then $x^T{M}_{e}x=0$ forces $x=0$, and an injective square map in finite dimensions is invertible; the same holds for $M_b$. Thus the system is exactly $z^{\mathrm{\prime }}=Az+f(t)$.

Let $g(s):=e^{-(s-t_0)A}f(s)$ and $H(t):={\int }_{t_0}^{t}g(s)ds$, interpreting each coordinate as a Riemann integral and using the oriented convention ${\int }_a^b=-{\int }_b^a$ when $b<a$ and ${\int }_a^a=0$. If $n_e+n_b=0$, then $g(s)$ and $H(t)$ are the unique empty vector for all $s,t$, so $H^{\mathrm{\prime }}=g$ vacuously. Otherwise, on every nondegenerate compact segment in $I$, each coordinate of $g$ is continuous, so [F3] gives integrability and [F4] defines $H$. For each coordinate $k$, the cut-additivity in [F6], together with the oriented convention, gives $H_k(t+h)-H_k(t)={\int }_t^{t+h}{g}_k(s)ds$ whenever $t,t+h\in I$: positive-oriented additivity follows by splitting at an interior point, and reversing endpoints or using a zero-length interval gives the other endpoint orders. Thus, for $h\mathrm{\ne }0$, $\frac{H_k(t+h)-H_k(t)}{h}-g_k(t)=\frac{1}{h}{\int }_t^{t+h}(g_k(s)-g_k(t))ds$, where [F7] evaluates the integral of the constant $g_k(t)$. Given $\epsilon >0$, continuity of $g_k$ at $t$ gives $\delta >0$ such that $\mathrm{\mid }{g}_k(s)-g_k(t)\mathrm{\mid }<\epsilon$ whenever $\mathrm{\mid }s-t\mathrm{\mid }<\delta$; since $I$ is open, reduce $\delta$ so $[t-\delta ,t+\delta ]\subset I$. If $0<\mathrm{\mid }h\mathrm{\mid }<\delta$, then on $J=[\min (t,t+h),\max (t,t+h)]$ the absolute-value estimate and monotonicity in [F6], followed by [F7], give $\mid \frac{1}{h}{\int }_t^{t+h}(g_k(s)-g_k(t))ds\mid \le \frac{1}{\mathrm{\mid }h\mathrm{\mid }}{\int }_J\mathrm{\mid }{g}_k(s)-g_k(t)\mathrm{\mid }ds\le \epsilon$. Therefore each coordinate has derivative $g_k(t)$ and $H^{\mathrm{\prime }}(t)=g(t)$ throughout $I$.

Define $z(t):=e^{(t-t_0)A}(z_0+H(t))$, where $z_0=(e_0,b_0)$. By [F2] and step 1.2, the product rule gives $z^{\mathrm{\prime }}=Az+f(t)$, since $e^{(t-t_0)A}{e}^{-(t-t_0)A}=I$ by the commuting-product identity in [F2]; also $z(t_0)=z_0$. Thus a $C^1$ solution exists on all of $I$. If $n_e+n_b=0$, every state is the unique empty vector, so the solution is unique. Otherwise, if $\tilde{z}$ is another solution with the same initial value, then $\xi (t):=e^{-(t-t_0)A}(\tilde{z}(t)-z(t))$ has derivative zero by [F2] and the two differential equations. Each coordinate of $\xi$ is constant by [F5] and vanishes at $t_0$, so $\tilde{z}=z$.

Differentiate $W=\frac{1}{2}(e^T{M}_{e}e+b^T{M}_{b}b)$ and substitute the system equations. Symmetry of $M_e,M_b$ gives $\dot{W}=e^T(C^T{M}_{b}b-Se-j)+b^T{M}_b(-Ce)=-e^{T}Se-e^{T}j$, because the real scalar cross terms $e^T{C}^T{M}_{b}b$ and $b^T{M}_{b}Ce$ are equal and cancel.

The magnetic quantity satisfies $(Db{)}^{\mathrm{\prime }}=D\dot{b}=-DCe=0$ by [F1]; each coordinate is constant by [F5], and its value at $t_0$ is $Db(t_0)$.

Since $q=-G^T{M}_{e}e$, differentiation and the first system equation give $\dot{q}=-G^T(C^T{M}_{b}b-Se-j)=-(CG{)}^T{M}_{b}b+G^T(Se+j)=G^T(Se+j)$, using $CG=0$ from [F1]. This proves all three identities.

For a zero-dimensional mass space, the unique square matrix is the identity, its own inverse, and its root. In positive dimension, $M_e$ and $M_b$ are symmetric and therefore self-adjoint in the standard Euclidean inner product; by [F1] each has an orthonormal eigenbasis. If $v$ is a unit eigenvector of either matrix with eigenvalue $\lambda$, positive definiteness in [F6] gives $\lambda =v^{T}Mv>0$. Define the square root to act by $\sqrt{\lambda }$ on the full $\lambda$-eigenspace, and its inverse to act by $1\mathrm{/}\sqrt{\lambda }$ there. These definitions are independent of a basis choice within an eigenspace, satisfy $(M^{1\mathrm{/}2}{)}^2=M$ and $M^{-1\mathrm{/}2}{M}^{1\mathrm{/}2}=I$, and agree with the Statement's root convention.

Let $x_n=M_e^{1\mathrm{/}2}{e}_n$ and $y_{n-1\mathrm{/}2}=M_b^{1\mathrm{/}2}{b}_{n-1\mathrm{/}2}$. Multiplying the first update in [F7] by $M_b^{1\mathrm{/}2}$ gives $y_{n+1\mathrm{/}2}=y_{n-1\mathrm{/}2}-hK{x}_n$. Multiplying the second update by $M_e^{1\mathrm{/}2}$ and using the root identities from step 1.1 gives $x_{n+1}=x_n+h{K}^T{y}_{n+1\mathrm{/}2}$.

Write $x=x_n$, $y=y_{n-1\mathrm{/}2}$, $y^{\mathrm{\prime }}=y_{n+1\mathrm{/}2}$, and $x^{\mathrm{\prime }}=x_{n+1}$. From step 2.1, $x^{\mathrm{\prime }}=x+h{K}^T{y}^{\mathrm{\prime }}$ and $y^{\mathrm{\prime }}=y-hKx$. Expanding gives $\parallel x^{\mathrm{\prime }}{\parallel }_2^2=\parallel x{\parallel }_2^2+2h(y^{\mathrm{\prime }}{)}^{T}Kx+h^2\parallel K^T{y}^{\mathrm{\prime }}{\parallel }_2^2$, $\parallel y^{\mathrm{\prime }}{\parallel }_2^2=\parallel y{\parallel }_2^2-2h{y}^{T}Kx+h^2\parallel Kx{\parallel }_2^2$, and $-h(y^{\mathrm{\prime }}{)}^{T}K{x}^{\mathrm{\prime }}=-h(y^{\mathrm{\prime }}{)}^{T}Kx-h^2\parallel K^T{y}^{\mathrm{\prime }}{\parallel }_2^2$.

Suppose $h\parallel K{\parallel }_2>2$. By the zero-domain clause of [F3], the domain is nonzero; [F2]–[F4] then give unit vectors $v,u$ for $\sigma =\parallel K{\parallel }_2$, with $Kv=\sigma u$ and $K^{T}u=\sigma v$. The recurrence in step 2.1 sends states $x=av$, $y=bu$ to states of the same form, with coefficient matrix $T_{\sigma }:=\left(\begin{array}{cc}1-h^2{\sigma }^2& h\sigma \\ -h\sigma & 1\end{array}\right)$. Its trace is $2-h^2{\sigma }^2<-2$ and determinant is $1$, so its characteristic polynomial has two distinct real negative reciprocal roots, one with modulus greater than $1$. A nonzero eigenvector for that root therefore has its coefficient norm multiplied by a factor greater than $1$ at every update, proving the supercritical growth claim.

Adding the identities in step 3.1 cancels the $\parallel K^T{y}^{\mathrm{\prime }}{\parallel }_2^2$ terms; substituting $y^{\mathrm{\prime }}=y-hKx$ in the remainder gives $\parallel x^{\mathrm{\prime }}{\parallel }_2^2+\parallel y^{\mathrm{\prime }}{\parallel }_2^2-h(y^{\mathrm{\prime }}{)}^{T}K{x}^{\mathrm{\prime }}=\parallel x{\parallel }_2^2+\parallel y{\parallel }_2^2-h{y}^{T}Kx$. Therefore $E_h(x^{\mathrm{\prime }},y^{\mathrm{\prime }})=E_h(x,y)$ at each update.

By [F5] and [F3]–[F4], $\mathrm{\mid }{y}^{T}Kx\mathrm{\mid }\le \parallel y{\parallel }_2\parallel Kx{\parallel }_2\le \parallel K{\parallel }_2\parallel x{\parallel }_2\parallel y{\parallel }_2$. Since $2\parallel x{\parallel }_2\parallel y{\parallel }_2\le \parallel x{\parallel }_2^2+\parallel y{\parallel }_2^2$, the cross term in $E_h$ has absolute value at most $\frac{1}{2}h\parallel K{\parallel }_2(\parallel x{\parallel }_2^2+\parallel y{\parallel }_2^2)$. Substitution gives both displayed energy bounds, whose lower coefficient is strictly positive when $h\parallel K{\parallel }_2<2$.

For every update, [F7] and [F6] give $D{b}_{n+1\mathrm{/}2}=D{b}_{n-1\mathrm{/}2}-hDC{e}_n=D{b}_{n-1\mathrm{/}2}$, so induction preserves $Db$. The loop in [F7] has exactly the finite index set $0,\dots ,N-1$; each exact mass solve is unique by [F6] and step 1.1. It therefore terminates after $N$ updates, including the case $N=0$.