Global ODE and Log-density Integration Error

Put $a:=1+Lh$ and $b:=C{h}^{p+1}$. By induction on $n$, the recurrence gives $e_n\le a^n{e}_0+b{\sum }_{j=0}^{n-1}{a}^j$ for $0\le n\le N$: the sum is empty at $n=0$, where the bound is equality in the first term, and substituting the bound at $n$ into $e_{n+1}\le a{e}_n+b$ appends the term $b{a}^n$.

If $L>0$, then $a-1=Lh$ and [F6, F7] give $(a-1){\sum }_{j<n}{a}^j={\sum }_{j<n}(a^{j+1}-a^j)=a^n-1$, so ${\sum }_{j<n}{a}^j=(a^n-1)\mathrm{/}(Lh)$. By [F2], $a=1+Lh\le \exp (Lh)$. If $n=0$ then $a^n=\exp (Lh{)}^n=1$; if $n\ge 1$, [F3] gives $a^n\le \exp (Lh{)}^n$. Induction using [F1] and [F4] gives $\exp (Lh{)}^n=\exp (nLh)$ for every natural $n$. Since $n\le N$, $nLh\le LT$, and [F5] gives $a^n\le \exp (LT)$. Substitution in step 1.1, using [F8] to simplify the real-power factor, yields $e_n\le (C_0\exp (LT)+C\frac{\exp (LT)-1}{L})h^p$. If $L=0$, then $a=1$ and step 1.1 gives $e_n\le C_0{h}^p+nC{h}^{p+1}=C_0{h}^p+C{h}^p(nh)\le (C_0+CT)h^p$ by [F8], because $n\le N$ and $Nh=T$. Taking the maximum proves both cases.

For fixed $t$, the inverse-flow chain rule [F5] gives $D{\mathrm{\Phi }}_t^{-1}({\mathrm{\Phi }}_t(x_0))D{\mathrm{\Phi }}_t(x_0)=I$, so $D{\mathrm{\Phi }}_t(x_0)$ is invertible and [F9] gives $\det D{\mathrm{\Phi }}_t(x_0)\mathrm{\ne }0$. Thus the formula for $p_t$ is continuous and positive on $U_t$ and zero outside, hence Borel measurable on ${\mathbb{R}}^d$. For any Borel $E\subseteq {\mathbb{R}}^d$, apply [F3] to $g={\mathbf{1}}_{E\cap U_t}{p}_t{\mathrm{\mid }}_{U_t}$. Change of variables gives ${\int }_E{p}_t(y)dy={\int }_{E\cap U_t}{p}_t(y)dy={\int }_{{\mathrm{\Phi }}_t^{-1}(E\cap U_t)}{p}_0(x)dx={\mu }_0({\mathrm{\Phi }}_t^{-1}(E))={\mu }_t(E)$ because ${\mathrm{\Phi }}_t(U_0)=U_t$. Thus $p_t$ is a density of ${\mu }_t$ on all of ${\mathbb{R}}^d$.

Put $s:=T\mathrm{/}2$. For $z\in U_s$ define ${\mathrm{\Psi }}_t(z):={\mathrm{\Phi }}_t({\mathrm{\Phi }}_s^{-1}(z))$ for $0\le t\le T$. Each ${\mathrm{\Psi }}_t(z)$ is a solution of the same ODE on the common interval $[0,T]=[s-s,s+s]$, and ${\mathrm{\Psi }}_s(z)=z$. Apply [F4] to this family at the initial state $x(s)={\mathrm{\Phi }}_s(x_0)$ and initial time $s$: its derivative $V_t:=D_z{\mathrm{\Psi }}_t(x(s))$ satisfies $V_t^{\mathrm{\prime }}=A_t{V}_t$, where $A_t:=D_{x}f(t,{\mathrm{\Phi }}_t(x_0))$, and $V_s=I_d$, with one-sided derivatives at $0,T$. Since ${\mathrm{\Phi }}_t={\mathrm{\Psi }}_t\circ {\mathrm{\Phi }}_s$, the chain rule [F5] gives $J_t:=D{\mathrm{\Phi }}_t(x_0)=V_t{J}_s$, so $J_t^{\mathrm{\prime }}=A_t{J}_t$ on $[0,T]$ and $J_0=I_d$, again with one-sided endpoint derivatives. Thus $J_t$ is continuous. Separately, the inverse-flow chain rule [F5] and [F9] give $\det J_t\mathrm{\ne }0$ for each $t$; [F10] and continuity of $J_t$ give continuity of $t\mapsto \det J_t$. If this determinant were negative at any $t>0$, [F6] on $[0,t]$ with target $0$ would give a zero, contradicting nonsingularity; since $\det J_0=1$ by [F11], it is positive throughout. The density formula in step 1.1 therefore becomes $p_t({\mathrm{\Phi }}_t(x_0))=p_0(x_0)\mathrm{/}\det J_t$.

Fix an interior time, write $A:=D_{x}f(t,{\mathrm{\Phi }}_t(x_0))$ and $C:=J_t^{-1}{J}_t^{\mathrm{\prime }}$, and denote the columns of $J_t$ by $v_1,\dots ,v_d$. By the Leibniz formula [F8], differentiating each $d$-factor monomial with repeated binary product rules (induction on the number of factors) and distributing the finite sum gives $(\det J_t{)}^{\mathrm{\prime }}={\sum }_{\sigma \in S_d}\mathrm{sgn}(\sigma ){\sum }_{j=1}^d(v_j^{\mathrm{\prime }}{)}_{\sigma (j)}{\prod }_{i\mathrm{\ne }j}(v_i{)}_{\sigma (i)}={\sum }_{j=1}^d\det [v_1,\dots ,v_{j-1},v_j^{\mathrm{\prime }},v_{j+1},\dots ,v_d]$. Since $J_t^{\mathrm{\prime }}=J_{t}C$, if $C_j$ is column $j$ of $C$ then $v_j^{\mathrm{\prime }}={\sum }_{k=1}^d{C}_{kj}{v}_k$. Multilinearity [F11] expands the $j$th determinant into ${\sum }_k{C}_{kj}\det [v_1,\dots ,v_{j-1},v_k,v_{j+1},\dots ,v_d]$; for $k\mathrm{\ne }j$ the repeated columns make the determinant zero by alternation [F11], while $k=j$ gives $C_{jj}\det J_t$. Hence $(\det J_t{)}^{\mathrm{\prime }}=\det J_t\mathrm{tr}(C)$ and [F7] gives $(\log \det J_t{)}^{\mathrm{\prime }}=\mathrm{tr}(C)$. Since $C=J_t^{-1}A{J}_t$, $\mathrm{tr}(C)={\sum }_{i,j,k}(J_t^{-1}{)}_{ij}{A}_{jk}(J_t{)}_{ki}={\sum }_{j,k}{A}_{jk}(J_t{J}_t^{-1}{)}_{kj}={\sum }_j{A}_{jj}=\mathrm{tr}(A)$. The density formula in step 2.1 has constant numerator $p_0(x_0)$, so [F7] and the chain rule [F5] give $(\log p_t({\mathrm{\Phi }}_t(x_0)){)}^{\mathrm{\prime }}=-(\det J_t{)}^{\mathrm{\prime }}\mathrm{/}\det J_t=-\mathrm{tr}(D_{x}f(t,{\mathrm{\Phi }}_t(x_0)))$. This proves the identity for $0<t<T$.

At the endpoints, step 2.1 gives the variational equation for $J_t$ on $[0,T]$, and [F1] makes its endpoint equations one-sided: $J_{0+}^{\mathrm{\prime }}=A_0{J}_0$ and $J_{T-}^{\mathrm{\prime }}=A_T{J}_T$, where $A_r:=D_{x}f(r,{\mathrm{\Phi }}_r(x_0))$. To differentiate the determinant there, use its finite Leibniz polynomial [F8]: as $t\downarrow 0$, each matrix entry has expansion $J_{ij}(t)=J_{ij}(0)+t(J_{0+}^{\mathrm{\prime }}{)}_{ij}+o(t)$, and as $t\uparrow T$ it has expansion $J_{ij}(t)=J_{ij}(T)+(t-T)(J_{T-}^{\mathrm{\prime }}{)}_{ij}+o(\mathrm{\mid }t-T\mathrm{\mid })$. Expanding each finite product in that polynomial, the terms linear in the increment are exactly the binary product-rule terms and all terms with at least two increments are little-oh of the increment. Thus the same finite expansion as in step 3.1 gives $(\det J{)}_{0+}^{\mathrm{\prime }}=\det J_0\mathrm{tr}(J_0^{-1}{A}_0{J}_0)=\det J_0\mathrm{tr}(A_0)$ and $(\det J{)}_{T-}^{\mathrm{\prime }}=\det J_T\mathrm{tr}(J_T^{-1}{A}_T{J}_T)=\det J_T\mathrm{tr}(A_T)$. The determinants are positive by step 2.1, so one-sided differentiation of their logarithms [F7] gives $(\log \det J{)}_{0+}^{\mathrm{\prime }}=\mathrm{tr}(A_0)$ and $(\log \det J{)}_{T-}^{\mathrm{\prime }}=\mathrm{tr}(A_T)$. Using the density identity in step 2.1, $\log p_t({\mathrm{\Phi }}_t(x_0))=\log p_0(x_0)-\log \det J_t$, therefore its right derivative at $0$ and left derivative at $T$ are respectively $-\mathrm{tr}(A_0)$ and $-\mathrm{tr}(A_T)$. These are exactly the stated endpoint identities.

By [F1], each $y(t_n)$ is a defined reference state of the exact initial-value problem; [F2] identifies the two uniform inequalities in [A1] as the one-step stability and local-defect bounds. Put $e_n:=\mathrm{\parallel }{y}_n-y(t_n){\mathrm{\parallel }}_2$. The Euclidean triangle inequality [F3] and those bounds give $e_{n+1}\le \mathrm{\parallel }{\mathrm{\Psi }}_{h,n}(y_n)-{\mathrm{\Psi }}_{h,n}(y(t_n)){\mathrm{\parallel }}_2+\mathrm{\parallel }{\mathrm{\Psi }}_{h,n}(y(t_n))-y(t_{n+1}){\mathrm{\parallel }}_2\le (1+Lh)e_n+C{h}^{p+1}$.

Apply [L1] to the nonnegative sequence $(e_n{)}_{n=0}^N$ using the recurrence from step 1.1 and the initial bound in [A1]. The two cases $L>0$ and $L=0$ give exactly the displayed conclusions, with exponential notation as in [F4].

Let $e_n:=\mathrm{\mid }{\widehat{\mathrm{\ell }}}_n-{\mathrm{\ell }}_n\mathrm{\mid }$. By [F1], $\mathrm{\parallel }{x}_n-x(t_n){\mathrm{\parallel }}_2\le K_x{h}^p$. The exact and numerical scalar updates, the triangle inequality, and [A1] give $e_{n+1}\le e_n+\mathrm{\mid }{Q}_h(t_n,x_n)-Q_h(t_n,x(t_n))\mathrm{\mid }+\mathrm{\mid }{Q}_h(t_n,x(t_n))-I_n\mathrm{\mid }\le e_n+h{K}_q{K}_x{h}^p+D{h}^{m+1}$.

Iterating the recurrence from step 1.1 for $n$ steps yields $e_n\le e_0+n(h{K}_q{K}_x{h}^p+D{h}^{m+1})$. Since $n\le N=T\mathrm{/}h$ and $e_0\le C_{\mathrm{\ell },0}{h}^r$, this is at most $C_{\mathrm{\ell },0}{h}^r+T{K}_q{K}_x{h}^p+DT{h}^m$. Taking the maximum proves the claim.

By [F1], the exact log density along the trajectory solves the scalar IVP ${\mathrm{\ell }}^{\mathrm{\prime }}(t)=q(t,x(t))$ with $\mathrm{\ell }(0)=\log p_0(x(0))$. Applying [F3] on each step shows that its increment from $t_n$ to $t_{n+1}$ is exactly $I_n={\int }_{t_n}^{t_{n+1}}q(s,x(s))ds$, so the exact scalar grid sequence in [F2] is ${\mathrm{\ell }}_n=\log p_{t_n}(x(t_n))$.

Apply [F2] to the numerical state and scalar updates. The exact grid sequence is the log density by step 1.1, so the stated bound is precisely the numerical log-density error.

For fixed $t_n$, the open convex fiber $U_{t_n}$ contains $K$. Applying [F3] there gives $\mathrm{\parallel }f(t_n,u)-f(t_n,v){\mathrm{\parallel }}_2\le L_f\mathrm{\parallel }u-v{\mathrm{\parallel }}_2$ for $u,v\in K$. Hence the Euler state map satisfies $\mathrm{\parallel }u+hf(t_n,u)-v-hf(t_n,v){\mathrm{\parallel }}_2\le (1+L_{f}h)\mathrm{\parallel }u-v{\mathrm{\parallel }}_2$.

Let $r(s):=f(t_n+s,x(t_n+s))$ for $0\le s\le h$. By [F2], $r^{\mathrm{\prime }}(s)={\mathrm{\partial }}_{t}f(t_n+s,x(t_n+s))+D_{x}f(t_n+s,x(t_n+s))f(t_n+s,x(t_n+s))$, whose norm is at most $M_f+L_f{B}_f$. By [F4], $\mathrm{\parallel }r(s)-r(0){\mathrm{\parallel }}_2\le (M_f+L_f{B}_f)s$. Using [F1] and the integral norm bound in [F4], the exact-start Euler defect is ${\parallel x(t_{n+1})-x(t_n)-hf(t_n,x(t_n))\parallel }_2={\parallel {\int }_0^h(r(s)-r(0))ds\parallel }_2\le C_x{h}^2$.

Since $f\in C^2(U)$, $q=-{\mathrm{div}}_{x}f$ is $C^1$ on $U$. The open convex fiber $U_{t_n}$ contains $K$, so [F3] gives $\mathrm{\mid }hq(t_n,u)-hq(t_n,v)\mathrm{\mid }\le h{L}_q\mathrm{\parallel }u-v{\mathrm{\parallel }}_2$ for $u,v\in K$. Along the exact trajectory, [F2] gives $\frac{d}{ds}q(t_n+s,x(t_n+s))={\mathrm{\partial }}_{t}q+D_{x}qf$, with absolute value at most $M_q+L_q{B}_f$, hence $\mathrm{\mid }q(t_n+s,x(t_n+s))-q(t_n,x(t_n))\mathrm{\mid }\le (M_q+L_q{B}_f)s$ by [F4]. The left-Euler exact-start defect is $\mid {\int }_0^h(q(t_n+s,x(t_n+s))-q(t_n,x(t_n)))ds\mid$, which [F4] bounds by $(M_q+L_q{B}_f)h^2\mathrm{/}2=D_q{h}^2$.