In Kendall Atkinson's An Introduction to Numerical Analysis, the author proves that the explicit Euler method for solving an inhomogeneous IVP is stable with respect to perturbations to the vector field and the initial iterate. As an exercise, I have tried to formulate and prove an analogous result for the implicit Trapezoidal method. I'm just looking to make sure my theorem statement and proof are correct.
Theorem. Let $[x_0,b]$ be a finite interval, let $h > 0$, let $N \equiv N(h) := \lfloor \frac{b-x_0}{h} \rfloor$, and let $f: [x_0,b] \times \mathbb{R} \to \mathbb{R}$ be uniformly Lipschitz in its second argument with constant $K \geq 0$ (that is, $|f(x,y_1) - f(x,y_2)| \leq K|y_1 - y_2|$ for all $x \in [x_0,b]$ and all $y_1,y_2 \in \mathbb{R}$), and assume that $hK \leq 1$. Let $y_0 \in \mathbb{R}$ be a chosen initial iterate and let $y_n \equiv y_h(x_n), \; 1 \leq n \leq N$ be the iterates obtained for solving the IVP $Y'(x) = f(x,Y(x)), \, Y(x_0) = Y_0$ by the (implicit) Trapezoidal method, given by \begin{equation} \hspace{1.2cm} y_{n+1} = y_n + \frac{h}{2}\left[f(x_n,y_n) + f(x_{n+1},y_{n+1}) \right], \qquad n \geq 0 \hspace{1.2cm} (1) \end{equation}
Let $\delta:[x_0,b]^2 \to \mathbb{R}$ be a bounded continuous function (denoting a structural perturbation to the vector field $f$) and let $\epsilon \in \mathbb{R}$ be some initial perturbation. Define a sequence of perturbed iterates $(z_n)_{n \geq 0}$ by setting $z_0 := y_0 + \epsilon$ defining $z_n$ for $n \geq 1$ according to the recursive equation \begin{equation} \hspace{1cm} z_{n+1} = z_n + \frac{h}{2}\left[f(x_n,z_n) + f(x_{n+1},z_{n+1}) + \delta(x_n) \right]. \hspace{1cm} (2) \end{equation}
Then there exists constants $C_1,C_2$ independent of $h$ such that \begin{align*} \hspace{1.8cm} \max_{0 \leq n \leq N(h)} |z_n - y_n| \leq C_1 |\epsilon| + C_2 \| \delta \|_u. \hspace{1.8cm} (3) \end{align*}
(Here, $\|\delta\|_u := \sup\{|\delta(s_1,s_2)|: (s_1,s_2) \in [x_0,b]^2 \}.)$ Below is my proof attempt.
Proof: Let $d_n := |z_n - y_n|$ for all $n \geq 0$. Subtracting (1) from (2) gives \begin{align*} d_{n+1} = d_n + \frac{h}{2} \left[ f(x_n,z_n) - f(x_n,y_n) \right] + \frac{h}{2} \left[f(x_{n+1},z_{n+1}) - f(x_{n+1},y_{n+1}) \right] + \frac{h}{2} \delta(x_{n-1},x_n). \end{align*}
Then applying the Triangle Inequality and the Lipschitz condition on $f$ gives \begin{align*} d_{n+1} &\leq d_n + \frac{h}{2}K|z_n - y_n| + \frac{h}{2}K|z_{n+1}-y_{n+1}| + \frac{h}{2} \|\delta\|_u \\[3pt] &= (1 + \tfrac{hK}{2})d_n + \tfrac{hK}{2}d_{n+1} + \tfrac{h}{2}\|\delta\|_u, \end{align*}
which then implies \begin{align*} d_{n+1} \leq \left(\frac{1 + \frac{hK}{2}}{1-\frac{hK}{2}}\right) d_n + \frac{h}{2}\|\delta\|_u. \end{align*} (The sign on the inequality does not flip because $hK \leq 1$). Now using the fact that $\frac{1+x}{1-x} \leq 1 + 4x$ for $0 \leq x < 1/2$, we obtain the bound \begin{align*} d_{n+1} \leq (1 + 2hK)d_n + \frac{h}{2} \|\delta\|_u. \end{align*}
Now observe that \begin{align*} d_1 &\leq (1 + 2hK) |\epsilon| + \tfrac{h}{2} \|\delta\|_u \\[3pt] d_2 &\leq (1 + 2hK) \left[ (1 + 2hK) |\epsilon| + \tfrac{h}{2} \|\delta\|_u \right] + \tfrac{h}{2} \|\delta\|_u \\[3pt] &= (1 + 2hK)^2 |\epsilon| + \left[1 + (1 + 2hK) \right] \tfrac{h}{2} \|\delta\|_u \\[3pt] a_3 &\leq (1 + 2hK) \left[(1 + 2hK)^2 |\epsilon| + \left[1 + (1 + 2hK) \right] \tfrac{h}{2} \|\delta\|_u \right] + \tfrac{h}{2} \|\delta\|_u \\[3pt] &= (1 + 2hK)^3 |\epsilon| + \left[1 + (1 + 2hK) + (1 + 2hK)^2 \right] \tfrac{h}{2} \|\delta\|_u \\[5pt] &\; \, \vdots \\ a_n &\leq (1 + 2hK)^n |\epsilon| + \left[(1 + (1 + 2hK) + \cdots + (1 + 2hK)^{n-1} \right] \tfrac{h}{2} \|\delta\|_u \\[5pt] &= (1 + 2hK)^n |\epsilon| + \left[\frac{(1 + 2hK)^n - 1}{(1+2hK)-1} \right] \tfrac{h}{2} \|\delta\|_u \\[5pt] &= (1 + 2hK)^n |\epsilon| + \left[\frac{(1 + 2hK)^n - 1}{4K} \right] \|\delta\|_u. \end{align*} where the second to last equality is by the formula for a finite geometric series. Now using the facts that $d_n \leq a_n$ and $(1 + x)^n \leq e^{nx}$ for all integers $n \geq 0$ and all $x \geq -1$, we have \begin{align*} (1 + 2hK)^n \leq e^{2Knh} = e^{2K(x_n-x_0)} \leq e^{2K(b-x_0)}. \end{align*}
Thus, \begin{align*} d_n \leq e^{2K(b-x_0)} |\epsilon| + \frac{e^{2K(b-x_0)} - 1}{4K} \|\delta\|_u. \end{align*}
And since the RHS of the above inequality is independent of $n$, it follows that \begin{align*} \max_{0 \leq n \leq N(h)} d_n \leq C_1 |\epsilon| + C_2 \|\delta\|_u \end{align*}
where $C_1 := e^{2K(b-x_0)}$ and $C_2 := \frac{e^{2K(b-x_0)} - 1}{4K}. \qquad \square$
Does the proof look okay?
Also, this is probably not as important but: when adding a vector field perturbation term $\delta$ to an implicit method of the form $y_{n+1} = G(y_{n-p},\ldots,y_{n+1})$, are the arguments of $\delta$ generally chosen to be $x_{n-p},\ldots,x_n$? Or do we also include $x_{n+1}$? (Of course, if we always assume that $\delta$ is bounded then this doesn't affect the proof, but I'm just curious.)