Question

Q. How to vary size of h in Runge Kutta Fehlberg Method?

Answer 1

Answer on Question #76733 – Math – Quantitative Methods

Question

How to vary size of $h$ in Runge Kutta Fehlberg Method?

Solution

To approximate the solution to the 1st order IVP:

$$y' = f(x, y), \qquad y(x_0) = y_0$$

we seek:

$$y_{n+1} = y_n + h \sum_{i=1}^{s} b_i k_i + O(h^{s+1})$$

The Runge-Kutta-Fehlberg method is a one-step method with the approximations calculated using the Runge-Kutta method of order 4 and 5. For this method each step requires the use of the following six values:

$$k_1 = h \cdot f(x_k, y_k)$$ $$k_2 = h \cdot f\left(x_k + \frac{1}{4} h, y_k + \frac{1}{4} k_1\right)$$ $$k_3 = h \cdot f\left(x_k + \frac{3}{8} h, y_k + \frac{3}{32} k_1 + \frac{9}{32} k_2\right)$$ $$k_4 = h \cdot f\left(x_k + \frac{12}{13} h, y_k + \frac{1932}{2197} k_1 - \frac{7200}{2197} k_2 + \frac{7296}{2197} k_3\right)$$ $$k_5 = h \cdot f\left(x_k + h, y_k + \frac{439}{216} k_1 - 8 k_2 + \frac{3680}{513} k_3 - \frac{845}{4104} k_4\right)$$ $$k_6 = h \cdot f\left(x_k + \frac{1}{2} h, y_k - \frac{8}{27} k_1 + 2 k_2 - \frac{3544}{2565} k_3 + \frac{1859}{4104} k_4 - \frac{11}{40} k_5\right)$$

Then we calculate the approximation to the solution with the help of the method of the fourth order:

$$y_{k+1}^4 = y_k + \frac{25}{216} k_1 + \frac{1408}{2565} k_3 + \frac{2197}{4101} k_4 - \frac{1}{5} k_5, \quad Error = O(h^4)$$

And the approximation to the solution with the help of the method of the 5th order:

$$y_{k+1}^5 = y_k + \frac{16}{135} k_1 + \frac{6656}{12825} k_3 + \frac{28561}{56430} k_4 - \frac{9}{50} k_5 + \frac{2}{55} k_6, \quad Error = O(h^5)$$

At each step, two different approximations for the solution are made and compared (Question #76732 - Math - Quantitative Methods):

$$\varepsilon_{n+1} = \frac{1}{h} \left| y_{k+1}^4 - y_{k+1}^5 \right|, \quad \delta < \left(\frac{\varepsilon}{\varepsilon_{n+1}}\right)^{1/4}$$

The optimal step size is $(\delta \cdot h)$. The Runge-Kutta-Felberg method consists in the fact that at each step of the method the accuracy of the function is determined by the difference in values between the results of the methods RK-4 and RK-5. If they differ by no more than $\varepsilon$- the required accuracy, then the value is considered an approximate value of the function at the point at the considered step. Otherwise, in the considered step, the new step value and the new value of the function are recomputed with the subsequent error estimate (1).

Let's consider stepping on the example (Question #76730 - Math – Quantitative Methods). In order for inequality (1) to be satisfied, it is necessary to multiply the optimal step by the number $\alpha \in (0,1)$. When calculating a new step, use formula:

$$\delta = \left(\frac {\varepsilon}{2 \cdot \varepsilon_ {n + 1}}\right) ^ {1 / 4} \cong 0. 8 4 \cdot \left(\frac {\varepsilon}{\varepsilon_ {n + 1}}\right) ^ {\frac {1}{4}} (2)$$

In addition, the solution can impose restrictions on the range of step change to obtain the required number of values of the function at a given interval.

If the two answers are in close agreement $(\varepsilon_{n + 1}\leq \varepsilon)$, then $\varepsilon_{n + 1}$ is the error of a method that is 5-th order accurate, then if we replace $h$ by $\delta \cdot h$, the error is multiplied by $\delta^5$. To calculate the new step and reduce the calculation time, it is used:

$$\left\{ \begin{array}{l} \delta = \alpha \cdot \left(\frac {\varepsilon}{\varepsilon_ {n + 1}}\right) ^ {\frac {1}{4}} i f \varepsilon_ {n + 1} > \varepsilon \\ \delta = \alpha \cdot \left(\frac {\varepsilon}{\varepsilon_ {n + 1}}\right) ^ {\frac {1}{5}} i f \varepsilon_ {n + 1} \leq \varepsilon \end{array} \right.$$

In addition, it is possible to reduce the computation time if we analyze the inequality $(\varepsilon_{n + 1}\leq \varepsilon)$:

$$\left\{ \begin{array}{c} \delta = \alpha \cdot \left(\frac {\varepsilon}{\varepsilon_ {n + 1}}\right) ^ {\frac {1}{4}} i f \varepsilon_ {n + 1} > \varepsilon \\ \delta = \alpha \cdot \left(\frac {\varepsilon}{\varepsilon_ {n + 1}}\right) ^ {\frac {1}{5}} i f \varepsilon / 2 < \varepsilon_ {n + 1} \leq \varepsilon \\ \delta = 2 \cdot h i f \varepsilon_ {n + 1} \leq \varepsilon / 2 \end{array} \right.$$

The solution by the classical method (2) requires 90000 steps:

fT(t) := 35 - 30·e-0.3·t

$Z\coloneqq$ Runge45(y0,D,Hmin,Hmax)

$i\coloneqq 0\dots$ rows(Z)-1 timei $\coloneqq Z_{i,0}$ T= fT(timei)

$\operatorname{rows}(Z) = 9 \times 10^4$

step $i\coloneqq Z_{i,2}$ AT $i\coloneqq Z_{i,1}$ errRKF45 $i\coloneqq Z_{i,3}$ calcCNT $i\coloneqq Z_{i,4}$

$\left|\mathrm{T} - \mathrm{AT}\right| = 1.02013767067$

$\mathrm{iC} := \sum_{i} \operatorname{calcCNT}_{i} = 9 \times 10^{4}$

if we analyze the inequality $(\varepsilon_{n + 1}\leq \varepsilon)$ , then a solution with required accuracy will be obtained at 1528 points with the calculation of the coefficients 3781 times:

$$fT(t) := 35 - 30 \cdot e^{-0.3 \cdot t} \quad Z := \operatorname{Runge} 45\left(y_0, D, Hmin, Hmax\right)$$ $$i := 0.. \operatorname{rows}(Z) - 1 \quad \underline{\text{time}_i} := Z_{i,0} \quad T_i := fT(\text{time}_i)$$ $$\operatorname{rows}(Z) = 1.528 \times 10^3 \quad \operatorname{step}_i := Z_{i,2} \quad AT_i := Z_{i,1} \quad \operatorname{errRKF45}_i := Z_{i,3} \quad \operatorname{calcCNT}_i := Z_{i,4}$$ $$|T - AT| = 0.12835905265 \quad iC := \sum_i \operatorname{calcCNT}_i = 3.781 \times 10^3$$

Answer provided by https://www.AssignmentExpert.com