Worked Error Diagnosis for Initial-Value Solvers

Diagnosing a Numerical Solution

A computed value for a differential equation is not automatically a reliable value. A method may use too few Taylor terms, a step size may be too large, or the numerical errors may grow even when the exact solution is well behaved. Diagnosis means checking each of these possibilities rather than treating the output as self-validating.

We will work with an initial-value problem (IVP), which specifies a differential equation and a starting value. The goal is to approximate the solution at later points, while keeping track of how the method and step size affect the result.

Definition: An IVP has the form \(y'=f(t,y)\), \(y(t_0)=y_0\). A one-step method uses the current approximation \((t_n,y_n)\) to calculate the next one, \(y_{n+1}\), at \(t_{n+1}=t_n+h\), where \(h\) is the step size.

Choosing a Taylor-Based Method

The Taylor method of order two follows by expanding the solution about \(t_n\) and retaining terms through \(h^2\). For a sufficiently differentiable solution,

y(t_n+h) = y(t_n)+h y'(t_n)+\frac{h^2}{2}y''(t_n)+O(h^3).

Since \(y'=f(t,y)\), the chain rule gives \(y''=f_t+f_y y'=f_t+f_y f\), with the partial derivatives evaluated at \((t,y)\). Thus the second-order Taylor update is

y_{n+1} = y_n+h f(t_n,y_n) +\frac{h^2}{2}\bigl(f_t(t_n,y_n)+f_y(t_n,y_n)f(t_n,y_n)\bigr).

This method is a reasonable candidate when the required derivatives are available and a second-order approximation is adequate. If the derivatives are inconvenient, another one-step method may be preferable. Whichever method is selected, the step size must also be checked: high formal order does not guarantee that a particular step is reliable.

Worked Example: Computing Two Taylor Steps

Consider \(y'=t+y\), \(y(0)=1\). Use the second-order Taylor method with \(h=0.1\).

Here \(f(t,y)=t+y\), so \(f_t=1\) and \(f_y=1\). Therefore \(y''=1+t+y\), and the update is

y_{n+1} =y_n+0.1(t_n+y_n)+0.005(1+t_n+y_n).

At the starting point, \(t_0=0\), \(y_0=1\). The derivative values are \(y'_0=0+1=1\) and \(y''_0=1+0+1=2\). Hence

y_1=1+0.1(1)+0.005(2)=1.110, \qquad t_1=0.1.

At the next point, \(y'_1=0.1+1.110=1.210\) and \(y''_1=1+0.1+1.110=2.210\). Thus

y_2=1.110+0.1(1.210)+0.005(2.210) =1.110+0.121+0.01105 =1.24205, \qquad t_2=0.2.

For a check, the exact solution is \(y(t)=2e^t-t-1\): differentiating gives \(y'=2e^t-1=t+y\), and \(y(0)=1\). Therefore \(y(0.2)=2e^{0.2}-1.2\approx1.242805516\). The computed value differs by about \(0.000755516\). At \(t=0.1\), the exact value is approximately \(1.110340372\), so the first-step error is about \(0.000340372\). The approximations are close in this example, but closeness should be assessed against the accuracy needed.

Local Error and Accumulated Error

Local truncation error describes the error made in one step when that step starts from the exact solution value. Global error is the difference between the computed value and the exact solution after multiple steps. Local errors do not simply disappear: later steps can reduce or amplify the errors already present.

A Taylor remainder gives a direct local error bound when a bound on a higher derivative is available. For a second-order Taylor step from the exact value \(y(t_n)\), if \(|y'''(t)|\le M\) throughout \([t_n,t_n+h]\), Taylor's theorem gives

\left| y(t_n+h)-\left(y(t_n)+h y'(t_n)+\frac{h^2}{2}y''(t_n)\right) \right| \le \frac{M h^3}{6}.

This is a local bound; to assess several steps, we must also account for how errors propagate. The following derivation makes that point explicit.

Write the exact one-step relation.
For an exact starting value, the second-order Taylor step differs from the true next value by a remainder \(r_n\), with \(|r_n|\le Mh^3/6\).

Track the error in the computed step.
For \(y'=-y\), the Taylor update is multiplication by \(q=1-h+h^2/2\). If \(E_n=y_n-y(t_n)\), then \(E_{n+1}=qE_n-r_n\).

Bound the accumulated error.
Taking absolute values gives \(|E_{n+1}|\le |q||E_n|+Mh^3/6\). This shows both the size of each new local error and the effect of propagating earlier errors.

The recurrence in the last step follows because the numerical update applied to \(y_n\) is \(q y_n\), while the exact next value is \(q y(t_n)+r_n\). Subtracting gives \(E_{n+1}=qE_n-r_n\). This argument is complete for \(y'=-y\); for a different equation, the propagation factor must be assessed from that method's update rather than copied automatically.

Worked Example: Bounding Error for a Decaying Solution

Consider \(y'=-y\), \(y(0)=1\), with \(h=0.2\). Since \(y''=y\), the second-order update is \(y_{n+1}=y_n(1-h+h^2/2)=0.82y_n\). Thus

y_1=0.82,\qquad y_2=0.82(0.82)=0.6724.

The exact solution is \(y(t)=e^{-t}\), which satisfies the IVP by differentiation and substitution. At \(t=0.2\), \(e^{-0.2}\approx0.818730753\), giving error \(|0.82-0.818730753|\approx0.001269247\). At \(t=0.4\), \(e^{-0.4}\approx0.670320046\), giving error \(|0.6724-0.670320046|\approx0.002079954\).

For the local remainder, \(y'''=-e^{-t}\), so \(|y'''|\le1\) on these steps. Each local error is therefore at most \(0.2^3/6=0.001333333\). Here \(q=0.82\), so the two-step global error satisfies

|E_1|\le0.001333333,\qquad |E_2|\le0.82(0.001333333)+0.001333333 \approx0.002426667.

The actual error \(0.002079954\) is below this bound. The bound is useful because it is guaranteed under its stated derivative condition, though it need not be tight.

Estimating Error by Comparing Step Sizes

A derivative bound may be difficult to obtain. A practical alternative is to compare one step of size \(h\) with two steps of size \(h/2\), using the same method. If the method's leading global error behaves like a constant times \(h^p\), the difference between the coarse and fine results is often approximately \((2^p-1)\) times the fine result's error. Dividing by \(2^p-1\) then estimates that fine error. This is an estimate, not a guaranteed bound: it is most informative when the step sizes are already in the method's asymptotic error regime.

Worked Example: A Step-Doubling Error Estimate

For \(y'=-y\), \(y(0)=1\), compare one second-order Taylor step of length \(0.2\) with two steps of length \(0.1\). The coarse result is \(0.82\), as calculated above. With \(h=0.1\), the update multiplier is \(1-0.1+0.1^2/2=0.905\), so the fine result is

y_{\mathrm{fine}}=0.905^2=0.819025.

The difference is \(|0.82-0.819025|=0.000975\). For an order-two method, divide by \(2^2-1=3\) to estimate the error in the fine result:

\text{estimated fine error}=\frac{0.000975}{3}=0.000325.

The exact value is \(e^{-0.2}\approx0.818730753\), so the actual fine-result error is \(|0.819025-0.818730753|\approx0.000294247\). The estimate is close, but it is not identical; it should not be reported as a rigorous bound.

Diagnosing a Stability Problem

For the test equation \(y'=\lambda y\), an order-two Taylor step multiplies the current approximation by \(R(z)=1+z+z^2/2\), where \(z=\lambda h\). When \(\lambda<0\), the exact solution decays. A step with \(|R(z)|>1\), however, amplifies numerical values from one step to the next. This is a warning that the step size is unsuitable for that problem and method.

Worked Example: A Step That Turns Decay into Growth

Take \(y'=-10y\), \(y(0)=1\), and use the second-order Taylor method with \(h=0.3\). Here \(\lambda=-10\), so \(z=-3\) and

R(-3)=1-3+\frac{(-3)^2}{2}=2.5.

The successive values are \(y_1=2.5(1)=2.5\) and \(y_2=2.5(2.5)=6.25\). Yet the exact solution \(y(t)=e^{-10t}\) has \(y(0.3)=e^{-3}\approx0.0497871\) and \(y(0.6)=e^{-6}\approx0.00247875\). The computed values grow while the exact solution rapidly decays. The diagnosis is not merely “some truncation error”: this step lies outside the method's stable range for this test equation.

Reducing the step to \(h=0.05\) gives \(z=-0.5\) and \(R(-0.5)=1-0.5+0.125=0.625\), whose magnitude is less than one. Over six steps to \(t=0.3\), the numerical value is \(0.625^6\approx0.059604645\), compared with \(e^{-3}\approx0.0497871\). The values now decay, although a visible approximation error remains. Stability is necessary for sensible behavior here, but it does not by itself guarantee the desired accuracy.

A Common Step-Size Mistake

A frequent error is to omit the square on the step size in the second-order term. For \(y'=-y\), \(y(0)=1\), and \(h=0.2\), one might incorrectly calculate \(y_1=1-0.2+\frac12=1.3\). This is wrong because Taylor's formula has \(\frac{h^2}{2}y''\), not \(\frac12 y''\). Since \(y''(0)=1\), the correct calculation is

y_1=1+(0.2)(-1)+\frac{(0.2)^2}{2}(1) =1-0.2+0.02 =0.82.

Checking the powers of \(h\) is especially important when changing step size: the first-order term scales with \(h\), while the second-order term scales with \(h^2\).

Putting the Diagnosis Together

When an approximation looks questionable, check the method formula and each substituted value first. Then compare the requested accuracy with a local or global error estimate, if available. If successive approximations behave unlike the expected solution, examine stability and try a smaller step. Finally, use an exact solution, when available, or a separate step-size comparison to check the result. These checks answer different questions: arithmetic checks catch implementation errors, error estimates assess accuracy, and stability checks detect harmful amplification.

Key takeaway: A trustworthy numerical solution requires more than a correctly applied formula. The method, step size, error estimate, and observed behavior must agree with the IVP and the accuracy required.

Check Your Understanding

Use the definitions and calculations in this tutorial to answer the following questions.

For \(y'=f(t,y)\), write the second-order Taylor update using \(f_t\) and \(f_y\).
For \(y'=-y\), \(y(0)=1\), and \(h=0.1\), compute one second-order Taylor step and compare it with \(e^{-0.1}\).
What information is needed to turn a Taylor remainder into a local error bound?
Why is the step-doubling estimate an estimate rather than a guaranteed error bound?
For \(y'=-8y\) with a second-order Taylor step and \(h=0.3\), calculate \(R(-2.4)\). Does its magnitude indicate amplification or decay?