View source for Numerical Analysis/Computing the order of numerical methods for ODE's

==<big>Numerical methods for ordinary differential equations order computation</big>==



===<big>Introduction</big>===

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 <math>O(h^{(n+1)})\,</math>, 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 [[w:Ordinary Differential Equations|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 [[w:Runge Kutta| 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. 



===<big>Example 1.  Determination of the parameters to establish a second order Runge Kutta method</big>===
Let a single step numerical method for solving ODE of the type

{{NumBlk|:|<math> y'(t)=f(t,y(t)), y(t_0)=y_0, \,</math>|{{EquationRef|1.1}}}}

be given by the following [2,4] recurrence formula 

<math> y_{n+1}=y_n+h(a_1k_1+a_2k_2) \,</math>

with

<math>k_1=f(t_n,y_n)\,</math>

<math>k_2=f(t_n+ph,y_n+hqk_1)\,</math>

Using Taylor series expansion about t<sub>n</sub>, the following can be obtained

{{NumBlk|:|<math> 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)\,</math>|{{EquationRef|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 <math> y(t_n)=y_n, y(t_n+h)=\overline{y_{n+1}}\,</math>. Therefore, both <math>y_n\,</math> and <math>\overline{y_{n+1}}\,</math> are considered exact values of <math>y(t)\,</math> at <math>t=t_n\,</math> and <math>t=t_n+h\,</math>, respectively.

Then, 

{{NumBlk|:|<math> \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)\,</math>|{{EquationRef|1.3}}}}


Using the differential equation (1.1) and replacing the derivatives of y, (1.3) becomes


{{NumBlk|:|<math> \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),\,</math>|{{EquationRef|1.4}}}}


where 

{{NumBlk|:|<math>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),\,</math>|{{EquationRef|1.5}}}}



or in the compact form (1.5) is <math>f'(t_n,y_n)=\left (f_t+f_yf) \right)_{t_n,y_n},\,</math>


and

{{NumBlk|:|<math>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}\,</math>|{{EquationRef|1.6}}}}


or in the compact form (1.6) is <math>f''(t_n,y_n)=\left ( f_{tt}+2f_{ty}f+f_tf_y+f_{yy}f^2+f_yf_yf \right )_{t=t_n,y=y_n}\,</math>


After substituting the expressions for f′ and f″ into the Taylor series expansion (1.4), the following is obtained

{{NumBlk|:|<math> \overline{y_{n+1}}=y_n+hf(t_n,y_n)+\frac{h^2}{2}(f_t+f_yf)_{(t_n,y_n)}+\frac{h^3}{6}\left ( f_{tt}+2f_{ty}f+f_tf_y+f_{yy}f^2+f_yf_yf \right )_{(t_n,y_n)}+O(h^4)\,</math>|{{EquationRef|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 <math> f(t_n+ph,y_n+hqk_1)\,</math> term about the <math>(t_n,y_n)\,</math> point, as follows.

{{NumBlk|:|<math>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^2f^2 \right )_{t=t_n,y=y_n}+O(h^3)\,</math>|{{EquationRef|1.8}}}}


or in the compact symbolic form (1.8) is:

{{NumBlk|:|<math>f(t_n+ph,y_n+hqk_1)=\left (f+f_t hp+f_yfhq+\frac {h^2}{2}(f_{tt}p^2+2f_{ty}qpf+f_{yy}q^2f^2) \right )_{t=t_n,y=y_n}+O(h^3)\,</math>|{{EquationRef|1.9}}}}


By substituting (1.9) into (1.7), the following is obtained

{{NumBlk|:|<math>y_{n+1}=y_n+\left(hf(a_1+a_2)+ha_2(hpf_t+hqf_yf)+\frac {h^3a_2}{2}(f_{tt} p^2+2f_tf_yqpf+f_{yy}q^2f^2)\right)_{t=t_n,y=y_n}+\underbrace{ha_2O_2(h^3)}_{O_3(h^4)}\,</math>|{{EquationRef|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.


{{NumBlk|:|<math>\overline{y_{n+1}}-y_{n+1}=\left((1-a_1-a_2)hf+(1/2-pa_2)h^2f_t+(1/2-qa_2)h^2f_yf+A(h)\right)_{t=t_n,y=y_n}\,</math>|{{EquationRef|1.11}}}}

where

{{NumBlk|:|<math>A(h)=\left( h^3(1/3(f_{tt}+2f_{ty}f+f_tf_y+f_{yy}f^2+f_yf_yf)-\frac{a_2}{2}(f_{tt}p^2+2f_tf_yqpf+f_{yy}q^2f^2))\right )_{t=t_n,y=y_n}+O_1(h^4)-O_3(h^4)\,</math>|{{EquationRef|1.12}}}}

It can be noticed that the multiplier expression of <math>h^3\,</math> 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 <math>O(h^3)\,</math>. In order to achieve the local truncation error of the third order, all terms in the error expression containing <math> h^0, h^1, h^2\,</math> must be eliminated, by adjusting the 4 parameters. By satisfying the equations (1.13-1.15), the local error is of <math>O(h^3)\,</math>, and, consequently, the global error is of <math>O(h^2)\,</math> order. Therefore, the resulting method is a second order method.


The equation set that must satisfied is the following 

{{NumBlk|:|<math> a_1+a_2=1, \,</math>|{{EquationRef|1.13}}}}

{{NumBlk|:|<math> pa_2=1/2, \,</math>|{{EquationRef|1.14}}}}

{{NumBlk|:|<math> qa_2=1/2. \,</math>|{{EquationRef|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 <math>a_1=a_2\,</math>, then <math>a_1=a_2=1/2, p=q=1\,</math>

The method equation (1.2) for these parameter becomes

<math> y_{n+1}=y_n+\frac{h}{2}(k_1+k_2) \,</math>

with

<math> k_1=f(t_n,y_n)\,</math>

<math>k_2=f(t_n+h,y_n+hk_1)\,</math>

or in a compact form

<math> y_{n+1}=y_n+\frac{h}{2}(f(t_n,y_n)+f(t_n+h,y_n+hf(t_n,y_n))) \,</math>


*(case 2) If <math> a_1=0\,</math> is chosen (notice that <math> a_2=0\,</math> cannot be selected, because of the last two equations), then


<math> y_{n+1}=y_n+hk_2 \,</math>


<math>k_1=f(t_n,y_n)\,</math>

<math>k_2=f(t_n+1/2h,y_n+1/2hk_1)\,</math>


or in a compact form


<math> y_{n+1}=y_n+hf(t_n+1/2h,y_n+1/2hf(t_n,y_n)) \,</math>

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


Let the recurrence equation of a method be given by the following of Runge Kutta type with three slope evaluations at each step


{{NumBlk|:|<math> y_{n+1}=y_n+h(a_1k_1+a_2k_2+a_3k_3) \,</math>|{{EquationRef|2.1}}}}


with


<math> k_1=f(t_n,y_n),\,</math>

<math> k_2=f(t_n+p_1h,y_n+q_{11}k_1),\,</math>

<math> k_3=f(t_n+p_2h,y_n+q_{21}k_1+q_{22}k_2),\,</math>


Taylor series expansion of <math>y(t_n+h)\,</math> about <math>t_n\,</math> 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.


{{NumBlk|:|<math> \overline{y_{n+1}}=y_n+hf(t_n,y_n)+\frac{h^2}{2}(f_t+f_yf)_{(t_n,y_n)}+\frac{h^3}{6}\left ( f_{tt}+2f_{ty}f+f_tf_y+f_{yy}f^2+f_yf_yf \right )_{(t_n,y_n)}+O(h^4)\,</math>|{{EquationRef|2.2}}}}


The Taylor expansion of the terms in (2.2) is shown up to <math>O(h^4)\,</math>, rather then up to <math>O(h^5)\,</math>, 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).


<div class="NavFrame"><div class="NavHead">
The Taylor expansion (equation 2.3)

</div><div class="NavContent" style="text-align:center">

<math>y_{n+1}=y_n+ha_1f(t_n,y_n)+ha_2(f+p_1hf_t+q_{11}hf_yf+\frac {h^2}{2}(f_{tt}p_1^2+f_{yy}q_{11}^2f^2+2f_{ty}p_1q_{11}f)+O_2(h^3))+\,</math>

<math>ha_3(f+p_2hf_t+f_y\left(q_{21}hf+q_{22}h(f+p_1hf_t+q_{11}hf_yf+O_3(h^2))\right)+\frac {1}{2}(f_{tt}(p_2h)^2+f_{yy}h^2(q_{21}f+\,</math>

<math>q_{22}(f+\underbrace{p_1hf_t+q_{11}hf_yf+O_4(h^2)}_{O_5(h)}))^2+2f_{ty}p_2h^2(q_{21}f+q_{22}(f+\underbrace{p_1hf_t+q_{11}hf_yf+O_4(h^2)}_{O_5(h)}))))\,</math>


</div></div><div class="NavEnd">&nbsp;</div>


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 <math>h^4\,</math>.

{{NumBlk|:|<math>y_{n+1}=y_n+h(a_1+a_2+a_3)f(t_n,y_n)+h^2(a_2(p_1f_t+q_{11}f_yf)+a_3(p_2f_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}^2f^2+2f_{ty}p_1q_{11}f)+\,</math>


<math>+a_3(q_{22}(p_1f_tf_y+q_{11}f_y^2f)+\frac{1}{2}(f_{tt}p_2^2+f_{yy}(q_{21}+q_{22})^2f^2+2f_{ty}p_2f(q_{21}+q_{22}))))_{(t_n,y_n)}+\,</math>


<math>+\underbrace{[hO_2(h^3)+h^2O_3(h^2)(f_yq_{22}a_3)+f_{yy}q_{22}a_3fh^3O_5(h)+a_3f_{ty}p_2q_{22}h^3O_5(h)]_{(t_n,y_n)}}_{O_6(h^4)}\,</math>|{{EquationRef|2,4}}}}


By comparing the two expressions (2,4) and (2,2), the following system of equations is obtained.

{{NumBlk|:|<math>a_1+a_2+a_3=1\,</math>|{{EquationRef|2.5}}}}

{{NumBlk|:|<math>p_1a_2+p_2a_3=1/2\,</math>|{{EquationRef|2.6}}}}

{{NumBlk|:|<math>q_{11}a_2+a_3(q_{21}+q_{22})=1/2\,</math>|{{EquationRef|2.7}}}}

{{NumBlk|:|<math>a_2p_1^2+a_1p_2^2=1/3\,</math>|{{EquationRef|2.8}}}}

{{NumBlk|:|<math>a_2p_1q_{11}+a_3p_2(q_{21}+q_{22})=1/3\,</math>|{{EquationRef|2.9}}}}

{{NumBlk|:|<math>a_3q_{22}p_1=1/6\,</math>|{{EquationRef|2.10}}}}

{{NumBlk|:|<math>a_2q_{11}^2+a_3(q_{21}+q_{22})^2=1/3\,</math>|{{EquationRef|2.11}}}}

{{NumBlk|:|<math>a_3q_{22}q_{11}=1/6\,</math>|{{EquationRef|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 <math>q_{11}=p_1\,</math>. Similarly, by substracting (2.11) from (2.9) equation, we see that <math>p_2=q_{21}+q_{22}\,</math>. 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 <math>p_2=1, q_{11}=1/2\,</math>, then we obtain the following recurrence equation.

{{NumBlk|:|<math>y_{n+1}=y_n+\frac{h}{6}(k_1+4k_2+k_3)\,</math>|{{EquationRef|2.13}}}}

where

<math>k_1=f(t_n,y_n)\,</math>

<math>k_2=f(t_n+1/2h,y_n+1/2hk_1)\,</math>

<math>k_3=f(t_n+h,y_n-hk_1+2hk_2)\,</math>



The recurrence equation (2.13) is known Runge Kutta third order method [1,3] ([[w:List_of_Runge%E2%80%93Kutta_methods|List of Runge–Kutta methods]]), which indicates that our approach was correct.


===<big>Quiz - Method order computations, local and global truncation error</big>===
<quiz display=simple>
{What is the main importance of having a higher ODE method order 
|type="()"}
-It always better to increase the accuracy of the solution
+It requires fewer steps for the same accuracy restriction
-We need the method to be better consistent
-We need the method to be more stable
-It provides better margin of stability

{The number of slope predictions (number of k's) in an RK method is related to the order of the method as follows. The order of the method is:
|type="()"}

-exactly equal to the number slope predictions
-equal to number of k's (slope evaluations) +1
+upper bounded by the number of k's
-not related to the number of k's

{An n-th order of the method means that
|type="()"}
-the global truncation error is of the order <math>O(h^{n+1})\,</math>
-the local truncation error is of the order <math>O(h^{n})\,</math>
+the global truncation error is of the order <math>O(h^{n})\,</math>
-the maximum of the global and local error is n-th order

{The order of the local truncation error is related to the global truncation error in the following way
|type="()"}
-they are about equal
-they are not related
-the global is one higher order than the local
+the local is one higher order than the global

{Is the following true? The global truncation error is always less then the local, since the errors partially cancel out over steps.
|type="()"}
-true
+false

{What is the highest order that we can achieve with a Runge Kutta type method by using 5 k's?
|type="()"}
-We can achieve 5-th order
-That method is not used, since it is not consistent.
+We can achieve 4-th order
-We can achieve even higher order, it is sixth order accurate.

{How is the order of single step methods related to consistency of the method?
|type="()"}
+At least first order of the method indicates that the method is consistent
-To be consistent, a method must be second order or higher 
-not related at all, a single step method can be inconsistent and be of any order.

{What is the main "tool" to prove the order of a numerical method for solving ODE's of type <math>y'=f(t,y(t))\,</math>?
|type="()"}
-Solving a specific "stiff" problem <math>y'=f(t,y(t))=-Ky\,</math> and showing that it converges for specific real positive K
-Taylor series expansion for one variable
+Taylor series expansion for one and two variables
-Solving the ODE equation analytically and then numerically

{If we show that all terms in a numerical method for ODE recurrence equation match the terms in the Taylor series expansion up to the terms with <math>h^n\,</math>, but we do not analyze whether further cancellation of the terms exist, the following is true
|type="()"}
-The method is of the order n
-The method is of the order n, since the global error is of the order <math>O(h^{n})\,</math>
+The order of the method could be higher than n-1, but it is at least n-1

{If we use an explicit n-th (n>1) order accurate multistep method with s steps (s>2) to solve an ODE, but we use an n-1 order method to calculate the missing initial points that we need to start using the multistep method, what is the global order of accuracy of our calculation?
|type="()"}
-There is no such combined method, we cannot do that
-The overall method is n-th order, since we just use the n-1 order method for few initial points
+The overall method is of the order n-1, since the initial error from the lower accurate method remains in the cummulative error (global)
-The order of the method could be either n or n-1
-The order of the method is something between n-1 and n, just like the RK 45 method

{What do we get as an order of error from the following expression

<math>\frac{1}{h}O(h^2)+10O(h^3)\,</math>?

|type="()"}

-Something between <math>O(h^2)\,</math> and <math>O(h^3)\,</math>
-<math>O(h^2)\,</math>
-<math>O(h^3)\,</math>
+<math>O(h)\,</math>

{What do we get as an order of error from the following expression

<math>sin(h^2)\frac{1}{h}O(h^2)\,</math>?

|type="()"}

-Something between <math>O(h^2)\,</math> and <math>O(h)\,</math>
-<math>O(h^2)\,</math>
+<math>O(h^3)\,</math>
-<math>O(h)\,</math>

{What do we get as an order of error from the following expression

<math>e^hcos(h)\frac{1}{h}O(h^2)\,</math>?

|type="()"}

-It must be of order <math>O(h^2)\,</math>
-<math>O(h^3)\,</math>
+<math>O(h)\,</math>

</quiz>


===<big>Conclusion</big>===

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 <math>y_{n+1}\,</math> using the previous value <math>y_n\,</math> and the right hand side function <math> f(t,y(t))\,</math> evaluated at <math>t_n\,</math>. 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 <math>t_n\,</math>. Hence, the function f appearing in the recurrence equation as an evaluation at time different than <math>t_n\,</math>, need to be replaced with a truncated Taylor series expansion about <math>(t_n,y_n)\,</math>.

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.


===References===

[1] A list of the Runge - Kutta methods at Wikipedia [[w:List_of_Runge%E2%80%93Kutta_methods|List of Runge–Kutta methods]], on November 5, 2010

[2] Ian Jacques and Colin Judd, "Numerical Analysis," Chapman and Hall, New York, 1987

[3] Topic Runge Kutta methods on Wikipedia [[w:Runge%E2%80%93Kutta_methods|Runge Kutta methods]], on November 5, 2010

[4] John H. Mathews and Kurtis K. Fink, "Numerical Methods Using Matlab," 4th Edition, Prentice-Hall Inc., 2004



[[Category:Numerical analysis]]
[[Category:Methods]]
[[Category:Examples]]
[[Category:Exercises]]
[[Category:Quizzes]]
[[Category:Ordinary differential equations]]