Question 10.12: Determining the step size for given accuracy. (a) Develop th......
Determining the step size for given accuracy.
(a) Develop the formula that can be used for calculating the step size for a required accuracy for the fourth-order Runge-Kutta method.
(b) For the given ODE, apply the formula from part (a) to find the magnitude of the step size such that the local truncation error will be less than 10^{-6}.
\frac{d y}{d x}=-\frac{y}{1+x^{2}}, \text { with } \quad y(0)=1, \text { for } 0<x<1 (10.189)
(a) Since the local truncation error for the fourth-order Runge-Kutta method is O\left(h^{5}\right), the truncation error in a solution with a step size of h can be expressed by (analogous to Eq. (10.182)):
e_{h}^{T R}=A h^{4} (10.182)
e_{h}^{T R}=A h^{5} (10.190)
When the step size is halved, the local truncation error when the independent variable advances two steps is:
2 e_{h / 2}^{T R}=\frac{A h^{5}}{16} (10.191)
Subtracting the last two equations gives a relation analogous to Eq. (10.184):
e_{h}^{T R}-2 e_{h / 2}^{T R}=A h^{4}-\frac{A h^{4}}{8}=\frac{7}{8} A h^{4} (10.184)
e_{h}^{T R}-2 e_{h / 2}^{T R}=A h^{5}-\frac{A h^{5}}{16}=\frac{15}{16} A h^{5} (10.192)
By using Eqs. (10.185) and (10.186),
\begin{aligned} &y_h^{N S}=y_h^{\text {Taylor }}\left(x_i+h\right)+e_h^{T R}& (10.185)\\ &y_{h / 2}^{N S}=y_h^{T a y l o r}\left(x_i+h\right)+2 e_{h / 2}^{T R} & (10.186) \end{aligned}
the difference in the truncation errors in Eq. (10.192) can be written in terms of the difference in the numerical solutions that are obtained with step sizes of h and h / 2 :
y_{h}^{N S}-y_{h/2}^{N S}=\frac{15}{16} A h^{5}, \quad \text { or } \quad A=\frac{16}{15} \frac{\left(y_{h}^{N S}-y_{h / 2}^{N S}\right)}{h^{5}} (10.193)
Once the two numerical solutions are obtained, the constant A in Eq. (10.193) can be determined. With A known, the step size for the required accuracy is calculated by using Eq. (10.190):
h=\left(\frac{\xi}{A}\right)^{1 / 5} (10.194)
where \xi is the desired error.
Script File
(b) To find the magnitude of the step size such that the local truncation error will be less than 10^{-6}, Eq. (10.189) is solved twice with the fourth-order Runge-Kutta method: once with h=1 and once with h=1 / 2. The solution of Eq. (10.189) is performed in the Command Window of MATLAB by using the user-defined function odeRK4 that was developed in Example 10-6:
>> [x, y] = odeRK4 (@Chap10Exmp22ODE, 0, 1, 1, 1)
From this solution y_{h}^{N S}=0.45666666666667.
>> [x, y] = odeRK4 (@Chap10Exmp22ODE, 0, 1, 0. 5, 1)
From this solution y_{h / 2}^{N S}=0.45599730642061.
The argument Chap10Exmp12ODE in the user-defined function odeRK4 is the name of the following user-defined function that calculates the value of d y / d x :
function dydx = Chap10Exmp12ODE (x, y)
dydx = -y/ (1 + x^2) ;
The results from the numerical solutions are used in Eq. (10.193) for determining the value of the constant A :
A=\frac{16}{15} \frac{(0.45666667-0.45599730)}{1^{5}}=7.14 \times 10^{-4}
Now that A is known, the magnitude of the step size such that the local truncation error will be less than 10^{-6} can be calculated with Eq. (10.194):
h=\left(\frac{10^{-6}}{7.14 \times 10^{-4}}\right)^{1 / 5}=0.2687
Thus, the step size must be smaller than h=0.2687 for solving Eq. (10.189) with the fourth-order Runge-Kutta method for the local truncation error to be smaller than 10^{-6}.
To check the result, Eq. (10.189) is solved again using h=0.2, which is smaller than the value obtained above.
>> [x,y]=odeRK4(@Chap10Exmp12ODE,0,1,0.2,1);
>> y(6)
ans =
0.45594036941322
This numerical solution gives y=0.45594036941322.
The exact solution of Eq. (10.189) is y=0.45593812776600.
The difference between the two solutions is 2.24 \times 10^{-6} which is larger than but on the order of 10^{-6}. However, the step size of h was calculated for a required local truncation error, and the last solution also contains a global truncation error (round-off errors are likely negligible since MATLAB uses double precision).