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

Question

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 \begin{equation} \hspace{3.7cm} Y'(x) = f(x,Y(x)), \quad Y(x_0) = Y_0. \hspace{3.8cm} (1) \end{equation}

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): \begin{equation} \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) \end{equation}

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 \begin{equation} \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) \end{equation}

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 \begin{equation} 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) \end{equation}

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:

Answer 1

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, \begin{equation} 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). \end{equation} 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, \begin{equation} \frac{1}{2}[Y^{(3)}(x_n) + Y^{(3)}(x_{n+1})] - Y^{(3)}(\eta) = O(h). \end{equation} 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

\begin{equation} 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). \end{equation}

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$

Answer 2

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.