# Derivation of an asymptotic error formula for the Trapezoidal method for IVPs

I am trying to prove the following theorem, which is a paraphrased version of Exercise 6.18 in Kendall Atkinson's An Introduction to Numerical Analysis:

Theorem. Let $$[x_0,b]$$ be a finite interval, let $$h > 0$$, and let $$N \equiv N(h) := \lfloor \frac{b-x_0}{h} \rfloor$$. Let $$f \in C^3([x_0,b] \times \mathbb{R},\mathbb{R})$$ be uniformly Lipschitz in its second argument with constant $$K \geq 0$$. Assume $$f_{yy}$$ is bounded and that $$hK \leq 1$$. Let $$Y:[x_0,b] \to \mathbb{R}$$ be the unique solution to the IVP $$$$\hspace{3.7cm} Y'(x) = f(x,Y(x)), \quad Y(x_0) = Y_0. \hspace{3.8cm} (1)$$$$

Let $$x_n := x_0 + nh$$ for $$n=0,1,\ldots,N(h)$$, let $$y_0 \in \mathbb{R}$$, and let $$(y_h(x_n))_{n \geq 1} = (y_n)_{n \geq 1}$$ be defined according to the implicit Trapezoidal method applied to (1): $$$$\hspace{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{2.1cm} (2)$$$$

Let $$e_n := Y(x_n) - y_n$$ be the error at the $$n$$-th iterate, let $$\delta_0 \in \mathbb{R}$$, and assume that the initial error satisfies $$e_0 = \delta_0 h^2 + O(h^3)$$. Then \begin{align*} \hspace{5cm} e_n = D(x_n)h^2 + O(h^3) \hspace{5cm} (3) \end{align*} for all $$n \geq 0$$, where $$D:[x_0,b] \to \mathbb{R}$$ is the solution of the IVP $$$$\hspace{1.95cm} D'(x) = f_y(x,Y(x))D(x) - \frac{1}{12}Y^{(3)}(x), \quad D(x_0) = \delta_0. \hspace{2.1cm} (4)$$$$

Some remarks:

• The Lipschitz condition on $$f$$ is not really relevant to this proof; it's just there to guarantee the convergence of the Trapezoidal method.
• Atkinson proves a similar asymptotic error formula for the forward Euler method (p.352-354), which I have included at the end of this post for reference.
• Assuming the proof of the above theorem is analogous, I think it should follow these steps:

Step 1. Derive a recurrence relation for the error $$e_n$$ in terms of $$h$$, $$f$$ and $$Y$$, using equation (2) and the truncation error of the trapezoidal method.

Step 2. Identify the so-called principal part of the error, call it $$g_n$$, in the recurrence formula for $$e_n$$. (The principal part of the error consists of the terms that make the dominant contribution to the error.) Set up a new recurrence formula of the form $$g_{n+1} := F(g_n)$$ for the accumulation of the principal part of the error.

Step 3. By making an appropriate change of variables from $$g_n$$ to $$\delta_n := h^{-2} g_n$$ (or something like that), obtain a formula that is the Trapezoidal method applied to the IVP in (4), with $$\delta_{n+1}$$ playing the role of $$y_{n+1}$$: \begin{align*} \hspace{5mm} \delta_{n+1} = \delta_n + \tfrac{h}{2} \big[ f_y(x_n,Y(x_n)) \delta_n - \tfrac{1}{12}Y^{(3)}(x_n) + f_y(x_{n+1},Y(x_{n+1})) \delta_{n+1} - \tfrac{1}{12}Y^{(3)}(x_{n+1}) \big] \hspace{5mm} (5) \end{align*}

Step 4. Prove that $$g_n = D(x_n)h^2 + O(h^3)$$ and $$e_n - g_n = O(h^m)$$ for some $$m \geq 3$$. The first equality follows easily from Step 3: Since the Trapezoidal method is convergent as $$h \to 0$$, we have \begin{align*} \delta_n - D(x_n) = O(h). \end{align*} Multiplying through by $$h^2$$ then gives $$g_n - D(x_n)h^2 = O(h^3).$$ I think I should be able to prove the second equality. We'll then have \begin{align*} e_n &= g_n + [e_n - g_n] \\[3pt] &= D(x_n)h^2 + O(h^3) + O(h^m) \\[3pt] &= D(x_n)h^2 + O(h^3) \end{align*} as desired.

My main difficulty is in Step 3. Here is what I have come up with so far.

Step 1: Firstly, since $$f \in C^3([x_0,b] \times \mathbb{R},\mathbb{R})$$, we have $$Y \in C^4([x_0,b],\mathbb{R})$$ (by this post). Then from the local truncation error formula of the Trapezoidal rule (equation (6.5.1) in Atkinson), \begin{align*} Y(x_{n+1}) = Y(x_n) + \frac{h}{2}\big[f(x_n,Y(x_n)) + f(x_{n+1},Y(x_{n+1})) \big] - \frac{h^3}{12} Y^{(3)}(\xi_n) \end{align*}

for some $$\xi_n \in [x_n,x_{n+1}]$$. By Taylor's Theorem applied to the last term, $$Y^{(3)}(x_n) = Y^{(3)}(x_n) + Y^{(4)}(\xi_n)(x_n - \xi_n)$$ and so we can write $$$$Y(x_{n+1}) = Y(x_n) + \frac{h}{2}\big[f(x_n,Y(x_n)) + f(x_{n+1},Y(x_{n+1})) \big] - \frac{h^3}{12} Y^{(3)}(x_n) + O(h^4). \qquad (6)$$$$

To ease notation, let $$Y_n := Y(x_n)$$. Then subtracting (6) from (2) gives \begin{align*} e_{n+1} &= e_n + \frac{h}{2}\left[f(x_n,Y_n) - f(x_n,y_n) + f(x_{n+1},Y_{n+1}) - f(x_{n+1},y_{n+1}) \right] \\[4pt] &\quad - \frac{h^3}{12} Y^{(3)}(x_n) + O(h^4). \qquad (7) \end{align*}

Now applying Taylor's theorem to the functions $$y \mapsto f(x_n,y)$$, we get \begin{align*} f(x_n,y_n) = f(x_n,Y_n) + f_y(x_n,Y_n)(y_n - Y_n) + f_{yy}(x_n,\zeta_n) \frac{(y_n - Y_n)^2}{2} \end{align*} for some $$\zeta_n$$ between $$y_n$$ and $$Y_n$$. Rearranging, we get \begin{align*} f(x_n,Y_n) - f(x_n,y_n) = e_n f_y(x_n,Y_n) - \tfrac{1}{2}f_{yy}(x_n,\zeta_n)e_n^2. \end{align*}

Similarly, \begin{align*} f(x_{n+1},Y_{n+1}) - f(x_n,y_n) = e_n f_y(x_{n+1},Y_{n+1}) - \tfrac{1}{2}f_{yy}(x_{n+1},\zeta_{n+1})e_n^2 \end{align*}

for some $$\zeta_{n+1}$$ between $$y_{n+1}$$ and $$Y_{n+1}$$. Substituting the above into (7) gives \begin{align*} e_{n+1} &= e_n + \frac{h}{2} \left[ e_n f_y(x_n,Y_n) - \tfrac{1}{2}f_{yy}(x_n,\zeta_n)e_n^2 + e_n f_y(x_{n+1},Y_{n+1}) - \tfrac{1}{2}f_{yy}(x_{n+1},\zeta_{n+1})e_n^2 \right] \\[5pt] & \quad -\frac{h^3}{12} Y^{(3)}(x_n) + O(h^4) \\[5pt] & \hspace{-5mm} = \left[1 + \frac{h}{2}f_y(x_n,Y_n) + \frac{h}{2}f_y(x_{n+1},Y_{n+1}) \right] e_n - \frac{h^3}{12} Y^{(3)}(x_n) - \frac{h}{4}e_n^2 \left[f_{yy}(x_n,\zeta_n) + f_{yy}(x_{n+1},\zeta_{n+1}) \right] + O(h^4) \end{align*}

Step 2. We should have $$e_n = O(h^2)$$ for the Trapezoidal method (by (6.5.18) in Atkinson), so we have $$\frac{h}{4}e_n^2 \left[f_{yy}(x_n,\zeta_n) + f_{yy}(x_{n+1},\zeta_{n+1}) \right] = O(h^5)$$. Dropping this term and the $$O(h^4)$$ terms, we define
\begin{align*} \hspace{2cm} g_{n+1} = \left[1 + \frac{h}{2}f_y(x_n,Y_n) + \frac{h}{2}f_y(x_{n+1},Y_{n+1}) \right] g_n - \frac{h^3}{12} Y^{(3)}(x_n), \qquad n \geq 0 \hspace{2cm} (8) \end{align*}

with $$g_0 := \delta_0 h^2$$, to represent the principal part of the error.

Step 3. Now define a new sequence $$\delta_n := g_n/h^2$$. Then write (8) in terms of $$\delta_n$$: \begin{align*} h^2 \delta_{n+1} = \left[1 + \frac{h}{2}f_y(x_n,Y_n) + \frac{h}{2}f_y(x_{n+1},Y_{n+1}) \right] h^2 \delta_n - \frac{h^3}{12} Y^{(3)}(x_n) \end{align*}

Cancelling out $$h^2$$ gives and rearranging leads to the equation \begin{align*} \delta_{n+1} = \delta_n + \frac{h}{2}\left[f_y(x_n,Y_n) \delta_n + f_y(x_{n+1},Y_{n+1}) \delta_n - \frac{1}{6} Y^{(3)}(x_n) \right] \end{align*}

This is the point where I got stuck. How can arrive at equation (5)? Or did I go wrong in a previous step?

Any help would be greatly appreciated!

For reference, here is Atkinson's proof of a similar theorem for Euler's method:

• You just say that the obtained iteration is the Euler method applied to the claimed differential equation. The difference is of order one, also as claimed, the exact DE solution gives a leading error term. Commented Jul 17, 2023 at 19:04
• @LutzLehmann: Thanks for your comment. Wouldn't Euler's method applied to (4) look like $\delta_{n+1} = \delta_n + h \left[f_y(x_n, Y(x_n)) \delta_n - \frac{1}{12} Y^{(3)}(x_n) \right]$ ? I'm not sure how to obtain that from what I have... Commented Jul 17, 2023 at 19:12
• You are right, the connection back is more by the same method, the implicit trapezoidal method. The general idea remains the same, look at $\frac{e_{n+1}-e_n}{h}$ as difference quotient, formally take the limit $h\to 0$, see that the resulting DE is regular enough so that the difference equation remains close to the differential equation. Commented Jul 17, 2023 at 19:19
• @LutzLehmann: I'm afraid I don't quite follow...Perhaps you could spell out the details in an answer? Commented Jul 17, 2023 at 19:27

You can find that because of symmetry you can even get $$\begin{multline} Y(x_{n+1}) = Y(x_n) + \frac{h}{2}\bigl[f(x_n,Y(x_n)) + f(x_{n+1},Y(x_{n+1})) \bigr]\\ - \frac{h^3}{24} \bigl[Y^{(3)}(x_n)+Y^{(3)}(x_{n+1})\bigr]+O(h^5) \end{multline}$$

Now insert the formula for the numerical approximation and compute the differences $$e_n=Y(x_n)-y_n$$ to get $$\begin{multline} e_{n+1} = e_n + \frac{h}{2}\bigl[f_y(x_n,Y(x_n))e_n + f_y(x_{n+1},Y(x_{n+1}))e_{n+1}+O(e_n^2,e_{n+1}^2) \bigr]\\ - \frac{h^3}{24} \bigl[Y^{(3)}(x_n)+Y^{(3)}(x_{n+1})\bigr]+O(h^5) \end{multline}$$ So if we work under the assumption of the result, then $$he_n^2=O(h^5)$$, so that the higher-order terms of the Taylor expansion do not have influence on the claimed result. One could at first also only assume that $$e_n=O(h)$$, to then bootstrap to $$e_n=O(h^2)$$.

Now compare this formula with the numerical method in question to detect that it is, in its main terms, the trapezoidal method for the differential equation $$e'(x)=f_y(x,Y(x))e(x)-\frac{h^2}{12}Y'''(x).$$ This is linear in $$e$$, the coefficients are quite non-singular, so that one can set $$e(x)=c(x)h^2$$, where $$c'(x)=f_y(x,Y(x))c(x)-\frac1{12}Y'''(x)$$ does not depend on $$h$$, the function $$c$$ is thus invariant under $$h$$. Considering that the order gap in all approximations used is $$2$$, the result for the numerical error is $$e_n=c(x_n)h^2+O(h^4).$$

To close out the circular argument one has to use a true induction based on inequalities. Especially the Taylor approximation, at least in the remainder term, gets replaced by mean-value inequalities, chiefly $$|f(x_n,Y(x_n))-f(x_n,y_n)|=\left|\int_0^1f_y(x_n,y_n+se_n)e_n\,ds\right|\le L_f|e_n|,\\ |Y'''(\xi)|\le M_3,$$ to get $$|e_{n+1}|\le |e_n|+\frac{L_fh}2(|e_n|+|e_{n+1}|)+\frac{h^3M_3}{12}$$ This, after reordering, is a recursive inequality of the form $$a_{n+1}\le qa_n+b$$ with $$a_n=|e_n|$$, $$q=\frac{2+Lh}{2-Lh}$$, $$b=\frac{Mh^3}{6(2-Lh)}$$ with solution $$q^{-n}a_{n}\le q^{-n+1}a_{n-1}+q^{-n}b\le...\le a_0+q^{-n}(q^{n-1}+...+q+1)b, \\~\\ a_n\le q^na_0+\frac{q^n-1}{q-1}b.$$ Inserting back this gives $$|e_n|\le \left(\frac{2+Lh}{2-Lh}\right)^n\left[|e_0|+\frac{M_3h^2}{12L}\right].$$

Note that $$\frac{2+x}{2-x}=\exp(\ln(1+x/2)-\ln(1-x/2))=\exp(x+x^3/12+x^5/80+...)=e^x+O(x^3)$$ so that $$\left(\frac{2+Lh}{2-Lh}\right)^n=e^{Lhn}+O(nh^3)=e^{L(x_n-x_0)}+O(h^2)$$ is mainly a constant relative to $$h$$, but depending on $$x_n$$.

Thus if also $$e_0=O(h^2)$$, for instance if $$y_0=Y(x_0)$$ as is usual, this proves that also $$e_n=O(h^2)$$. Thus the derivation of the differential equation for the leading term is justified.

• Thanks for your answer! Could you explain how you obtained the last two terms of the first equation? Specifically, why is $-\frac{h^3}{12} Y^{(3)}(\xi_n) = - \frac{h^3}{24}\left[ Y^{(3)}(x_n) + Y^{(3)}(x_{n+1}) \right] + O(h^5)$? Commented Jul 20, 2023 at 16:37
• After a bit of reading, I found the corrected trapezoidal rule in Atkinson (p.255): $\int_{a}^{b} f(x)\,dx \approx \frac{h}{2}[f(a) + f(b)] - \frac{h^2}{12}[f'(b) - f'(a)]$. This looks very similar to what you used in the first equation, but I'm still confused because we have $f'(b) - f'(a)$ in the last term, rather than $f'(b) + f'(a)$.... Commented Jul 20, 2023 at 18:28
• That is not correct in general. here it results from the time symmetry of the method. You can express everything in Taylor expansions at $x_{n+1/2}$ and then the remaining derivatives in the leading terms by symmetric expressions in $x_n,x_{n+1}$. // Apply the trapezoidal method again to get $f'(b)-f'(a)=\frac{h}2(f''(a)+f''(b)]+O(h^3)$. Commented Jul 20, 2023 at 21:34
• I still don't follow. Could you spell out the details? Commented Jul 21, 2023 at 14:03
• I made a new post to address this question: math.stackexchange.com/questions/4740731/… Commented Jul 22, 2023 at 18:59

Here's a proof in my own words. (Credit to @Lutz Lehmann for the help.)

We begin with an essential lemma. (See this post for a stronger result.)

Lemma. Let $$Y \in C^4([a,b],\mathbb{R})$$, let $$x_n < x_{n+1}$$ be in $$[a,b]$$, and let $$h := x_{n+1} - x_n$$. Then \begin{align*} Y(x_{n+1}) = Y(x_n) + \frac{h}{2}\left[Y'(x_n) + Y'(x_{n+1})\right] - \frac{h^3}{24}\left[Y^{(3)}(x_n) + Y^{(3)}(x_{n+1}) \right] + O(h^4). \qquad (*) \end{align*}
Proof. The truncation error of the Trapezoidal method at $$x_n$$ is $$-\frac{h^3}{12} Y^{(3)}(\eta)$$, where $$\eta \in [x_n,x_{n+1}]$$. (See p.252 of Atkinson.) Thus, $$$$Y(x_{n+1}) = Y(x_n) + \frac{h}{2}\left[Y'(x_n) + Y'(x_{n+1})\right] - \frac{h^3}{12} Y^{(3)}(\eta).$$$$ Comparing the above with $$(*)$$, we see that it suffices to show \begin{align*} -\frac{h^3}{12}Y^{(3)}(\eta) = -\frac{h^3}{24}\left[Y^{(3)}(x_n) + Y^{(3)}(x_{n+1}) \right] + O(h^4), \end{align*} or equivalently, $$$$\frac{1}{2}[Y^{(3)}(x_n) + Y^{(3)}(x_{n+1})] - Y^{(3)}(\eta) = O(h).$$$$ To prove the above, note that $$\frac{1}{2}[Y^{(3)}(x_n) + Y^{(3)}(x_{n+1})]$$ lies between $$Y^{(3)}(x_n)$$ and $$Y^{(3)}(x_{n+1})$$. Since $$Y^{(3)}$$ is continuous, the Intermediate Value Theorem guarantees the existence of some $$c \in [x_n,x_{n+1}]$$ such that \begin{align*} \frac{1}{2}[Y^{(3)}(x_n) + Y^{(3)}(x_{n+1})] = Y^{(3)}(c). \end{align*}

And since $$Y^{(3)}$$ is continuously differentiable, Taylor's theorem then gives
\begin{align*} Y^{(3)}(c) = Y^{(3)}(\eta) + (c-\xi) Y^{(4)}(\xi) \end{align*} for some $$\xi$$ between $$c$$ and $$\eta$$. Since $$c$$ and $$\eta$$ are both in $$[x_n,x_{n+1}]$$, so is $$\xi$$. Therefore, $$|c - \xi| \leq h$$. Hence, $$Y^{(3)}(c) = Y^{(3)}(\eta) + O(h)$$, and so \begin{align*} \frac{1}{2}[Y^{(3)}(x_n) + Y^{(3)}(x_{n+1})] - Y^{(3)}(\eta) = Y^{(3)}(c) - Y^{(3)}(\eta) = O(h), \end{align*} completing the proof.

Proof of the Theorem. Firstly, since $$f \in C^3([x_0,b] \times \mathbb{R},\mathbb{R})$$, we have $$Y \in C^4([x_0,b],\mathbb{R})$$ (see this post). Then by the lemma and the fact that $$Y'(x) = f(x,Y(x))$$, we have

$$$$Y(x_{n+1}) = Y(x_n) + \frac{h}{2}\left[f(x_n,Y(x_n)) + f(x_{n+1},Y(x_{n+1})) \right] - \frac{h^3}{24}\left[Y^{(3)}(x_n) + Y^{(3)}(x_{n+1}) \right] + O(h^4).$$$$

Now let $$e_n := Y(x_n) - y_n$$ for all $$n$$. Subtracting the above equation from (2) yields the recursive formula \begin{align*} e_{n+1} &= e_n + \frac{h}{2}\left[f(x_n,Y(x_n)) - f(x_n,y_n) \right] + \frac{h}{2}\left[f(x_{n+1},Y(x_{n+1})) - f(x_{n+1},y_{n+1}) \right] \\[5pt] & \quad - \frac{h^3}{24}\left[Y^{(3)}(x_n) + Y^{(3)}(x_{n+1}) \right] + O(h^4). \end{align*}

Now by Taylor's theorem applied to the map $$y \mapsto f(x_n,y)$$ (expanded about $$Y(x_n)$$) we have \begin{align*} f(x_n,y_n) &= f(x_n,Y(x_n)) + f_y(x_n,Y(x_n))(y_n - Y_n) + f_{yy}(x_n,\zeta_n) \frac{(y_n - Y(x_n))^2}{2} \end{align*} for some $$\zeta_n$$ between $$y_n$$ and $$Y(x_n)$$. Rearranging, we get

\begin{align*} f(x_n,Y(x_n)) - f(x_n,y_n) = f_y(x_n,Y(x_n))e_n + \tfrac{1}{2}f_{yy}(x_n,\zeta_n)e_n^2. \end{align*}

Now plugging in the above formula into the recursive equation for $$e_n$$, we obtain \begin{align*} e_{n+1} &= e_n + \tfrac{h}{2}\left[f_y(x_n,Y(x_n))e_n + \tfrac{1}{2}f_{yy}(x_n,\zeta_n)e_n^2 \right] + \tfrac{h}{2}\left[f_y(x_{n+1},Y(x_{n+1}))e_{n+1} + \tfrac{1}{2}f_{yy}(x_{n+1},\zeta_{n+1})e_{n+1}^2 \right] \\[5pt] &\quad - \tfrac{h^3}{24}\left[Y^{(3)}(x_n) + Y^{(3)}(x_{n+1}) \right] + O(h^4) \\[5pt] &= e_n + \tfrac{h}{2}\left[f_y(x_n,Y(x_n))e_n - \tfrac{h^2}{12}Y^{(3)}(x_n) + f_y(x_{n+1},Y(x_{n+1}))e_{n+1} - \tfrac{h^2}{12} Y^{(3)}(x_{n+1}) \right] \\[5pt] &\quad + O(he_{n}^2) + O(he_{n+1}^2) + O(h^4). \end{align*}

According to Atkinson (p.370), the Trapezoidal errors satisfy \begin{align*} \hspace{2cm} \max_{0 \leq n \leq N(h)} |e_n| \leq e^{2K(b-x_0)}|e_0| + \left[\frac{e^{2K(b-x_0)} - 1}{K} \right] \left[\frac{h^2}{12} \|Y^{(3)} \|_{\infty} \right], \hspace{2cm} (*) \end{align*}

where $$\|Y^{(3)} \|_{\infty} := \max_{x \in [a,b]} |Y^{(3)}(x)|$$. This is a direct consequence of Theorem 6.6 (p.360-361) of Atkinson. Now since $$e_0 = O(h^2)$$ by assumption, the above inequality implies $$e_n = O(h^2)$$ for all $$n \geq 0$$. Therefore, $$he_n^2 = O(h^5)$$ for all $$n$$, and so

\begin{align*} e_{n+1} &= e_n + \frac{h}{2}\left[f_y(x_n,Y(x_n))e_n - \frac{h^2}{12}Y^{(3)}(x_n) + f_y(x_{n+1},Y(x_{n+1}))e_{n+1} - \frac{h^2}{12} Y^{(3)}(x_{n+1}) \right] + O(h^4). \end{align*}

Now define a new sequence $$(g_n)$$ representing the principal part of the error: Let $$g_0 := \delta_0 h^2$$ and set \begin{align*} g_{n+1} &= g_n + \frac{h}{2}\left[f_y(x_n,Y(x_n))g_n - \frac{h^2}{12}Y^{(3)}(x_n) + f_y(x_{n+1},Y(x_{n+1}))g_{n+1} - \frac{h^2}{12} Y^{(3)}(x_{n+1}) \right] \end{align*}

for all $$n \geq 0$$. Then define a sequence $$(\delta_n)_{n\geq 0}$$ by setting $$\delta_{n} = h^{-2} g_n$$ for all $$n \geq 0$$. Then we have $$g_n = h^2 \delta_n$$, and so the previous equation becomes \begin{align*} h^2 \delta_{n+1} &= h^2 \delta_n + \frac{h}{2}\left[f_y(x_n,Y(x_n))h^2 \delta_n - \frac{h^2}{12}Y^{(3)}(x_n) + f_y(x_{n+1},Y(x_{n+1})) h^2 \delta_{n+1} - \frac{h^2}{12} Y^{(3)}(x_{n+1}) \right]. \end{align*}

After cancelling out $$h^2$$ from both sides, we obtain \begin{align*} \delta_{n+1} &= \delta_n + \frac{h}{2}\left[f_y(x_n,Y(x_n)) \delta_n - \frac{1}{12}Y^{(3)}(x_n) + f_y(x_{n+1},Y(x_{n+1})) \delta_{n+1} - \frac{1}{12} Y^{(3)}(x_{n+1}) \right]. \end{align*}

This is precisely the Trapezoidal method applied to the IVP (4), with $$\delta_n$$ playing the role of $$y_n$$. Since the initial error is zero in this case (because $$D(x_0) = \delta_0$$), the inequality $$(*)$$ gives the following error estimate: \begin{align*} D(x_n) - \delta_n = O(h^2). \end{align*} Multiplying this by $$h^2$$ and rearranging gives \begin{align*} g_n = D(x_n)h^2 + O(h^4). \end{align*}

Now let $$k_n := e_n - g_n$$ for all $$n \geq 0$$. By using the recursive formulas for $$e_n$$ and $$g_n$$, we obtain a recursive formula for $$k_n$$: \begin{align*} k_{n+1} = k_n + \frac{h}{2}\left[f_{yy}(x_n,\zeta_n)e_n^2 + f_{yy}(x_{n+1},\zeta_{n+1})e_{n+1}^2 \right] + O(h^4). \end{align*}

Now let $$M := \sup\{|f_{yy}(x,y)|:(x,y) \in [a,b] \times \mathbb{R} \}$$. (We have $$M < \infty$$ because $$f_{yy}$$ is assumed to be bounded.) Then applying the Triangle Inequality, \begin{align*} |k_{n+1}| \leq |k_n| + \frac{Mh}{2}\left[e_n^2 + e_{n+1}^2 \right]. \end{align*}

It follows from $$(*)$$ that $$\max_{0 \leq n \leq N(h)} |e_n^2| = O(h^4)$$; hence, there exists a constant $$C \geq 0$$ independent of $$h$$ and $$n$$ such that $$e_n^2 + e_{n+1}^2 \leq Ch^4$$ for all sufficiently small $$h$$. Therefore, \begin{align*} |k_{n+1}| \leq |k_n| + \tilde{C}h^5, \end{align*}

where $$\tilde{C} := \frac{MC}{2}.$$ It follows by induction that
\begin{align*} |k_n| &\leq |k_0| + n\tilde{C}h^5 \\[3pt] &\leq |k_0| + N(h) \tilde{C}h^5 \\[3pt] &\leq |k_0| + (b-x_0)\tilde{C}h^4 \\[3pt] &= |k_0| + O(h^4). \end{align*}

Now $$k_0 = e_0 - g_0 = [\delta_0 h^2 + O(h^3)] - \delta_0 h^2 = O(h^3)$$, and hence \begin{align*} k_n = O(h^3). \end{align*}

It follows that \begin{align*} e_n = g_n + k_n = \left[D(x_n)h^2 + O(h^4) \right] + O(h^3) = D(x_n)h^2 + O(h^3) \end{align*} as claimed. $$\qquad \square$$