Section 19.3 Estimating the accuracy of numeric ODE methods
Evaluating the accuracy of ODE solution methods is a more complex issue because they involve approximate computations based on the results of previous other approximate computations. A convenient model problem to use is \(y'=y\) with \(y(0)=1\text{.}\) Its exact solution at \(x= h\) is
The approximate solution \(y_1\) at \(x= h\) will be some other expression. The difference \(|e^h - y_1|\) gives us the local error of the method (for one step). The cumulative error can be estimated by multiplying the local error by the necessary number of steps to cover an interval \([a, b]\text{,}\) which is \((b-a)/h\text{.}\)
Euler's method has \(y_1 = 1 + h\text{,}\) hence \(|e^h - y_1| = h^2/2 + \dots\text{,}\) the local error is of second order. The cumulative error is of order \(1\text{.}\) This is consistent with what we know about left endpoint method of integration.
The above does not capture the entire picture. Although errors might accumulate, they do not always do that. When solving an equation like \(y' = -ry\) (exponential decay) we get \(y_k = (1-rh)^k y_0\text{.}\) As long as \(h \lt 2/r\text{,}\) we have \(|1-rh| \lt 1\) and so the approximate solution decays exponentially as it should. But if \(h \gt 2/r\text{,}\) the solution grows in magnitude instead of decreasing. This leads to the subject of stability of numerical ODE methods, which is beyond our scope.