# Numerical Analysis/Computing the order of numerical methods for ODE's
## Numerical methods for ordinary differential equations order computation
### Introduction

The subject of this analysis is the order of accuracy of numerical methods for solving ordinary differential equations. When we mention "order of accuracy" for a particular numerical method, we usually mean the order of the global truncation error. The global error is the cumulative error in the numerical solution that is produced on an interval that we need the ODE to solve on. A local truncation error is the error in the numerical solution (the round off errors produced by computing are excluded) generated at a particular step, when the previous step solution is considered exact (although it is not exact, unless the previous step is the boundary or initial condition). The local truncation error is usually denoted as $O(h^{(n+1)})\,$, where n is the order of accuracy on the entire interval of the differential equation and h is the step size (for fixed step size discretization). The importance of the order of a numerical method is reflected in the fact that we need fewer steps to numerically solve an ODE with a prescribed error tolerance, if we use a higher order numerical method. This issue becomes very important when we try to solve an ODE that describes a relatively fast changing physical proccess, such as an internal combustion in the IC machines, vibrations in structures, etc.
The following two examples show how to determine the order of a numerical method. The first example is about a set of Runge Kutta methods of the second order. It shows how we can derive the conditions on the parameters that we introduce in the general form of the method and the possible choices for the parameters. The second example can be used as an exercise, since some steps are hidden and can be seen by selecting the marked button to show the missing steps that complete the derivation. The example shows how the conditions are derived for a Runge Kutta type general form with slope evaluations at three different locations within the current step. The main discussion is about the comparison of the Taylor series expansion and the corresponding numerical method recurrence equation. Based on this comparison, term by term, we conclude on the order of the local truncation error, which is then used to make a statement about the order of the method. The examples are followed by a concept quiz, which is convenient for a user to review the theoretical details that are used in the examples. Finally, we end up the analysis with a conclusion.

### Example 1. Determination of the parameters to establish a second order Runge Kutta method

Let a single step numerical method for solving ODE of the type
    $y'(t)=f(t,y(t)),y(t_{0})=y_{0},\,$         (1.1)
be given by the following [2,4] recurrence formula
$y_{n+1}=y_{n}+h(a_{1}k_{1}+a_{2}k_{2})\,$
with
$k_{1}=f(t_{n},y_{n})\,$
$k_{2}=f(t_{n}+ph,y_{n}+hqk_{1})\,$
Using Taylor series expansion about t$_n$, the following can be obtained
    $y(t_{n}+h)=y(t_{n})+y'(t_{n})h+y''(t_{n})h^{2}/2+y'''(t_{n})h^{3}/3!+O(h^{4})\,$         (1.2)
Assuming that the previous point in (1.2) is exact (which we may because we are analyzing the local truncation error, i.e. the error generated at the last step), we will denote $y(t_{n})=y_{n},y(t_{n}+h)={\overline {y_{n+1}}}\,$. Therefore, both $y_{n}\,$ and ${\overline {y_{n+1}}}\,$ are considered exact values of $y(t)\,$ at $t=t_{n}\,$ and $t=t_{n}+h\,$, respectively.
Then,
    ${\overline {y_{n+1}}}=y_{n}+y'(t_{n})h+y''(t_{n})h^{2}/2+y'''(t_{n})h^{3}/3!+O(h^{4})\,$         (1.3)
Using the differential equation (1.1) and replacing the derivatives of y, (1.3) becomes
    ${\overline {y_{n+1}}}=y_{n}+f(t_{n},y_{n})h+f'(t_{n},y_{n})h^{2}/2+f''(t_{n},y_{n})h^{3}/3!+O(h^{4}),\,$         (1.4)
where
    $f'(t_{n},y_{n})={\frac {\partial f}{\partial t}}(t_{n},y_{n})+{\frac {\partial f}{\partial y}}(t_{n},y_{n})f(t_{n},y_{n}),\,$         (1.5)
or in the compact form (1.5) is $f'(t_{n},y_{n})=\left(f_{t}+f_{y}f)\right)_{t_{n},y_{n}},\,$
and
    $f''(t_{n},y_{n})=\left({\frac {\partial ^{2}f}{\partial t^{2}}}+{\frac {\partial ^{2}f}{\partial y\partial t}}f+{\frac {\partial ^{2}f}{\partial y\partial t}}f+{\frac {\partial f}{\partial y}}{\frac {\partial f}{\partial t}}+{\frac {\partial ^{2}f}{\partial y^{2}}}f^{2}+{\frac {\partial f}{\partial y}}{\frac {\partial f}{\partial y}}f\right)_{t=t_{n},y=y_{n}}\,$         (1.6)
or in the compact form (1.6) is $f''(t_{n},y_{n})=\left(f_{tt}+2f_{ty}f+f_{t}f_{y}+f_{yy}f^{2}+f_{y}f_{y}f\right)_{t=t_{n},y=y_{n}}\,$
After substituting the expressions for f′ and f″ into the Taylor series expansion (1.4), the following is obtained
    ${\overline {y_{n+1}}}=y_{n}+hf(t_{n},y_{n})+{\frac {h^{2}}{2}}(f_{t}+f_{y}f)_{(t_{n},y_{n})}+{\frac {h^{3}}{6}}\left(f_{tt}+2f_{ty}f+f_{t}f_{y}+f_{yy}f^{2}+f_{y}f_{y}f\right)_{(t_{n},y_{n})}+O(h^{4})\,$         (1.7)
Now, our objective is to adjust the method's recurrence equation (1.2), such that it can be compared, term by term, with the equation obtained using the Taylor series expansion (1.7). This step requires the Taylor series expansion of the $f(t_{n}+ph,y_{n}+hqk_{1})\,$ term about the $(t_{n},y_{n})\,$ point, as follows.
    $f(t_{n}+ph,y_{n}+hqk_{1})=f(t_{n},y_{n})+{\frac {\partial f}{\partial t}}(t_{n},y_{n})(hp)+{\frac {\partial f}{\partial y}}(t_{n},y_{n})(qhf(t_{n},y_{n}))+{\frac {h^{2}}{2}}\left({\frac {\partial ^{2}f}{\partial t^{2}}}p^{2}+2{\frac {\partial ^{2}f}{\partial y\partial t}}qpf+{\frac {\partial ^{2}f}{\partial y^{2}}}q^{2}f^{2}\right)_{t=t_{n},y=y_{n}}+O(h^{3})\,$         (1.8)
or in the compact symbolic form (1.8) is:
    $f(t_{n}+ph,y_{n}+hqk_{1})=\left(f+f_{t}hp+f_{y}fhq+{\frac {h^{2}}{2}}(f_{tt}p^{2}+2f_{ty}qpf+f_{yy}q^{2}f^{2})\right)_{t=t_{n},y=y_{n}}+O(h^{3})\,$         (1.9)
By substituting (1.9) into (1.7), the following is obtained
    $y_{n+1}=y_{n}+\left(hf(a_{1}+a_{2})+ha_{2}(hpf_{t}+hqf_{y}f)+{\frac {h^{3}a_{2}}{2}}(f_{tt}p^{2}+2f_{t}f_{y}qpf+f_{yy}q^{2}f^{2})\right)_{t=t_{n},y=y_{n}}+\underbrace {ha_{2}O_{2}(h^{3})} _{O_{3}(h^{4})}\,$         (1.10)
Finally, the method's equation (1.10) and the exact one step solution (1.7) can be compared. By substracting (1.10) from (10.7), the following error expression is obtained.
    ${\overline {y_{n+1}}}-y_{n+1}=\left((1-a_{1}-a_{2})hf+(1/2-pa_{2})h^{2}f_{t}+(1/2-qa_{2})h^{2}f_{y}f+A(h)\right)_{t=t_{n},y=y_{n}}\,$         (1.11)
where
    $A(h)=\left(h^{3}(1/3(f_{tt}+2f_{ty}f+f_{t}f_{y}+f_{yy}f^{2}+f_{y}f_{y}f)-{\frac {a_{2}}{2}}(f_{tt}p^{2}+2f_{t}f_{y}qpf+f_{yy}q^{2}f^{2}))\right)_{t=t_{n},y=y_{n}}+O_{1}(h^{4})-O_{3}(h^{4})\,$         (1.12)
It can be noticed that the multiplier expression of $h^{3}\,$ in (1.12) does not cancel out in (1.11) for any combination of the parameters, since there are terms that appear only once without any parameter. This means that the dominant term in the error expression cannot be of higher order than $O(h^{3})\,$. In order to achieve the local truncation error of the third order, all terms in the error expression containing $h^{0},h^{1},h^{2}\,$ must be eliminated, by adjusting the 4 parameters. By satisfying the equations (1.13-1.15), the local error is of $O(h^{3})\,$, and, consequently, the global error is of $O(h^{2})\,$ order. Therefore, the resulting method is a second order method.
The equation set that must satisfied is the following
    $a_{1}+a_{2}=1,\,$         (1.13)
    $pa_{2}=1/2,\,$         (1.14)
    $qa_{2}=1/2.\,$         (1.15)
There are three equations (1.13,1.14,1.15) but with four unknown parameters. This means that one parameter must be chosen. Now, it is tempting to ask whether the extra parameter (additional degree of freedom) could be used to cancel next term in the error expression. However, we have already seen that it is not possible in this case, due to the additional term that are not associated (multiplied) with a parameter.
Since one parameter can be chosen, then there is non-unique form of a second order method.
* (case 1) Choose $a_{1}=a_{2}\,$, then $a_{1}=a_{2}=1/2,p=q=1\,$
The method equation (1.2) for these parameter becomes
$y_{n+1}=y_{n}+{\frac {h}{2}}(k_{1}+k_{2})\,$
with
$k_{1}=f(t_{n},y_{n})\,$
$k_{2}=f(t_{n}+h,y_{n}+hk_{1})\,$
or in a compact form
$y_{n+1}=y_{n}+{\frac {h}{2}}(f(t_{n},y_{n})+f(t_{n}+h,y_{n}+hf(t_{n},y_{n})))\,$
* (case 2) If $a_{1}=0\,$ is chosen (notice that $a_{2}=0\,$ cannot be selected, because of the last two equations), then
$y_{n+1}=y_{n}+hk_{2}\,$
$k_{1}=f(t_{n},y_{n})\,$
$k_{2}=f(t_{n}+1/2h,y_{n}+1/2hk_{1})\,$
or in a compact form
$y_{n+1}=y_{n}+hf(t_{n}+1/2h,y_{n}+1/2hf(t_{n},y_{n}))\,$

### Example/Exercise 2. A single step ODE numerical method order computing with three slope evaluations (Runge Kutta 3-rd order)

Let the recurrence equation of a method be given by the following of Runge Kutta type with three slope evaluations at each step
    $y_{n+1}=y_{n}+h(a_{1}k_{1}+a_{2}k_{2}+a_{3}k_{3})\,$         (2.1)
with
$k_{1}=f(t_{n},y_{n}),\,$
$k_{2}=f(t_{n}+p_{1}h,y_{n}+q_{11}k_{1}),\,$
$k_{3}=f(t_{n}+p_{2}h,y_{n}+q_{21}k_{1}+q_{22}k_{2}),\,$
Taylor series expansion of $y(t_{n}+h)\,$ about $t_{n}\,$ is the same as in Example 1. Therefore, we will just use the final expression (1.7), since the procedure of the derivation is the same. For convenience, the final expression is repeated, which is going to be a reference equation for the comparison with the method's recurrence equation. Since the formulas for the given form of recurrence equation will get complicated, we will use the compact symbolic notation for the derivatives, which is shown in Example 1.
    ${\overline {y_{n+1}}}=y_{n}+hf(t_{n},y_{n})+{\frac {h^{2}}{2}}(f_{t}+f_{y}f)_{(t_{n},y_{n})}+{\frac {h^{3}}{6}}\left(f_{tt}+2f_{ty}f+f_{t}f_{y}+f_{yy}f^{2}+f_{y}f_{y}f\right)_{(t_{n},y_{n})}+O(h^{4})\,$         (2.2)
The Taylor expansion of the terms in (2.2) is shown up to $O(h^{4})\,$, rather then up to $O(h^{5})\,$, as we should in order to check that eventually next higher order terms cancel out, but we will assume that the method cannot achieve better local accuracy then fourth order, or equivalently, the global error of the third order. This will save us getting into the third level expansion of the two variable function f, which has 18 terms and would not be appropriate due to its length (even if the compact symbolic notation is used).
After we prepared the Taylor series expansion, we need to adjust the method's recurrence equation such that it can be compared with the Taylor series (2.2).
The Taylor expansion (equation 2.3)
$y_{n+1}=y_{n}+ha_{1}f(t_{n},y_{n})+ha_{2}(f+p_{1}hf_{t}+q_{11}hf_{y}f+{\frac {h^{2}}{2}}(f_{tt}p_{1}^{2}+f_{yy}q_{11}^{2}f^{2}+2f_{ty}p_{1}q_{11}f)+O_{2}(h^{3}))+\,$
$ha_{3}(f+p_{2}hf_{t}+f_{y}\left(q_{21}hf+q_{22}h(f+p_{1}hf_{t}+q_{11}hf_{y}f+O_{3}(h^{2}))\right)+{\frac {1}{2}}(f_{tt}(p_{2}h)^{2}+f_{yy}h^{2}(q_{21}f+\,$
$q_{22}(f+\underbrace {p_{1}hf_{t}+q_{11}hf_{y}f+O_{4}(h^{2})} _{O_{5}(h)}))^{2}+2f_{ty}p_{2}h^{2}(q_{21}f+q_{22}(f+\underbrace {p_{1}hf_{t}+q_{11}hf_{y}f+O_{4}(h^{2})} _{O_{5}(h)}))))\,$
Now, we need to group the terms in the similar way they are grouped in the Taylor series (2.2), such that we can establish the conditions on the parameters that will yield the same terms as in the Taylor expansion up to the terms containing $h^{4}\,$.
    $y_{n+1}=y_{n}+h(a_{1}+a_{2}+a_{3})f(t_{n},y_{n})+h^{2}(a_{2}(p_{1}f_{t}+q_{11}f_{y}f)+a_{3}(p_{2}f_{t}+q_{21}ff_{y}+q_{22}ff_{y}))_{(t_{n},y_{n})}+h^{3}({\frac {a_{2}}{2}}(f_{tt}p_{1}^{2}+f_{yy}q_{11}^{2}f^{2}+2f_{ty}p_{1}q_{11}f)+\,$ $+a_{3}(q_{22}(p_{1}f_{t}f_{y}+q_{11}f_{y}^{2}f)+{\frac {1}{2}}(f_{tt}p_{2}^{2}+f_{yy}(q_{21}+q_{22})^{2}f^{2}+2f_{ty}p_{2}f(q_{21}+q_{22}))))_{(t_{n},y_{n})}+\,$ $+\underbrace {[hO_{2}(h^{3})+h^{2}O_{3}(h^{2})(f_{y}q_{22}a_{3})+f_{yy}q_{22}a_{3}fh^{3}O_{5}(h)+a_{3}f_{ty}p_{2}q_{22}h^{3}O_{5}(h)]_{(t_{n},y_{n})}} _{O_{6}(h^{4})}\,$         (2,4)
By comparing the two expressions (2,4) and (2,2), the following system of equations is obtained.
    $a_{1}+a_{2}+a_{3}=1\,$         (2.5)
    $p_{1}a_{2}+p_{2}a_{3}=1/2\,$         (2.6)
    $q_{11}a_{2}+a_{3}(q_{21}+q_{22})=1/2\,$         (2.7)
    $a_{2}p_{1}^{2}+a_{1}p_{2}^{2}=1/3\,$         (2.8)
    $a_{2}p_{1}q_{11}+a_{3}p_{2}(q_{21}+q_{22})=1/3\,$         (2.9)
    $a_{3}q_{22}p_{1}=1/6\,$         (2.10)
    $a_{2}q_{11}^{2}+a_{3}(q_{21}+q_{22})^{2}=1/3\,$         (2.11)
    $a_{3}q_{22}q_{11}=1/6\,$         (2.12)
At the first glance, the system is closed, the number of equations is (2.5 through 2.12) which matches the number of undetermined parameters. However, only 6 equations are independent, the rest of them can be obtained from those 6 equations. By dividing (2.10) with (2.12), we can obtain that $q_{11}=p_{1}\,$. Similarly, by substracting (2.11) from (2.9) equation, we see that $p_{2}=q_{21}+q_{22}\,$. When we replace these two results into the rest of the equations, it is evident that the (2.6) and the (2.7) are the same, and (2.8) and the (2.9) equations are the same. Therefore, two equations can be obtained from other six, and we have to choose two variables in order to obtain a solution for the parameters.
For example, we can choose that $p_{2}=1,q_{11}=1/2\,$, then we obtain the following recurrence equation.
    $y_{n+1}=y_{n}+{\frac {h}{6}}(k_{1}+4k_{2}+k_{3})\,$         (2.13)
where
$k_{1}=f(t_{n},y_{n})\,$
$k_{2}=f(t_{n}+1/2h,y_{n}+1/2hk_{1})\,$
$k_{3}=f(t_{n}+h,y_{n}-hk_{1}+2hk_{2})\,$
The recurrence equation (2.13) is known Runge Kutta third order method [1,3] (List of Runge–Kutta methods), which indicates that our approach was correct.

### Conclusion

The main reason that we need to consider the order of a numerical method for solving ODE's is that it directly determines the step size i.e. the number of steps that we need to calculate the recurrence equation for. This directly influences the amount of computational work to obtain the prescribed accuracy.
The main approach for computing the order of a numerical method is that we first use Taylor series expansion for the exact differential equation to obtain the next value, namely $y_{n+1}\,$ using the previous value $y_{n}\,$ and the right hand side function $f(t,y(t))\,$ evaluated at $t_{n}\,$. We cannot say in advance up to which order we need to expand those terms in the Taylor series, since we are solving for the order of the method. The general rule is to expand the terms one order higher that we expect the method order is.
The second part of the method order computing is to write all terms in the given method recurrence equation in terms of the functions f and y evaluated at $t_{n}\,$. Hence, the function f appearing in the recurrence equation as an evaluation at time different than $t_{n}\,$, need to be replaced with a truncated Taylor series expansion about $(t_{n},y_{n})\,$.
The third part is to compare those two expressions obtained in the two forementioned ways. The order with respect to the powers of h, up to which the terms from the adjusted recurrence equation match the terms in the Taylor series expansion, represent the order of local truncation error. After we computed the order of the local truncation error, the global truncation error is one order less than the local error order.
Since the order of a numerical method for ODE is, by the definition, equal to the global error order, we only need to compute the local truncation error order and substract one from it.
There are some methods that do not have preciselly determined order. It is usually the case with the predictor corrector methods, like the Runge-Kutta method 45.