Numerical Analysis/stability of RK methods

Definitions

Let $\tau _{i}$ be the local truncation error, and k is the number of time steps, then:

1. The numerical method is consistent with a differential equation if

   $\lim _{h\to 0}\operatorname {max} |\tau _{i}|=0$ over $1\leqslant i\leqslant k$.

   according to this definition w:Euler's method is consistent.
2. A numerical method is said be convergent with respect to differential equation if

   $\lim _{h\to 0}|x(t_{i})-y_{i}|=0$ over $1\leqslant i\leqslant k$;

   where $y_{i}$ is the approximation for $x(t_{i})$.
3. A numerical method is stable if small change in the initial conditions or data, produce a correspondingly small change in the subsequent approximations.

Theorem: For an initial value problem

$x'=f(t,x)$ with $t\in [t_{0},t_{0}+\alpha ]$

and certain initial conditions, $(t_{0},x_{0})$, let us consider a numerical method of the form

$(y_{0}=x_{0})$ and $y_{i+1}=y_{i}+h\phi (t_{i},y_{i},h)$.

If there exists a value $h>0$ such that it is continuous on the iterative domain, $\Omega$ and if there exists an $L>0$ such that

$|\phi (t,y,h)-\phi (t,y^{*},h)|\leqslant L|y-y^{*}|$ for all $(t,y,h),\,(t,y^{*},h)\in \Omega$,

the method fulfills the w:Lipschitz condition, and it is stable and convergent if and only if it is consistent. That is,

$\phi (t,x,0)=f(t,x)$ for all $t\in \Omega$.

For a similar argument, one can deduce the following for multi- step methods:

1. The method is stable if and only if all roots, $\lambda$, of the characteristic polynomial satisfy

   $|\lambda |\leqslant 1$,

   and any root of

   $|\lambda |=1$

   is simple root.
2. one more result is that if the method is consistent with the differential equation, the method is stable if and only if it it is convergent.

see [1]

stability polynomials of Runge-Kutta methods

The w:Runge–Kutta methods are very useful in solving systems of differential equations, it has wide applications for the scientists and the engineers, as well as for the economical models, the recognized with their practical accuracy where we can use and get very good results and approximations when solving an ODE problem, RK has the general form :$y_{n+1}=y_{n}+h\sum _{i=1}^{s}b_{i}k_{i}$, where

$k_{1}=hf(t_{n},y_{n}),\,$

$k_{2}=hf(t_{n}+c_{2}h,y_{n}+a_{21}k_{1}),\,$

$k_{3}=hf(t_{n}+c_{3}h,y_{n}+a_{31}k_{1}+a_{32}k_{2}),\,$

$\vdots$

$k_{s}=hf(t_{n}+c_{s}h,y_{n}+a_{s1}k_{1}+a_{s2}k_{2}+\cdots +a_{s,s-1}k_{s-1}).$

such that

s is the number of stages,

a_ij (for 1 ≤ j < i ≤ s),

b_i (for i = 1, 2, ..., s)

and c_i (for i = 2, 3, ..., s).

Example:finding the stability polynomial for RK4's methods

for RK4" case$2_{a}$", which characterizes by $c_{2}={\frac {-1}{4}}$ which has the form: the stability region is found by applying the method to the linear test equation $y'=\lambda y$

$y_{n+1}=y_{n}+{\frac {1}{6}}k_{1}+0k_{2}+{\frac {2k_{3}}{3}}+{\frac {k_{4}}{6}}$

the stability region is found by applying the method to the linear test equation $y'=\lambda y$

$\displaystyle \ k_{1}=hf(t_{n},y_{n})$

$k_{2}=hf(t_{n}+{\frac {h}{2}},y_{n}-{\frac {k_{1}}{2}})$

$k_{3}=hf(t_{n}+{\frac {h}{2}},y_{n}+{\frac {3k_{1}}{4}}-{\frac {k_{2}}{4}})$

$k_{4}=hf(t_{n}+h,y_{n}-2k_{1}+k_{2}+2k_{3})$

using the linearized equation $f(t,y)=\lambda y$, and considering ${\hat {h}}=h\lambda$ , we get;

$k_{1}={\hat {h}}y_{n}$

$k_{2}={\hat {h}}(1-{\frac {\hat {h}}{2}})y_{n}$

$k_{3}={\hat {h}}(1+{\frac {3{\hat {h}}}{4}}-{\frac {\hat {h}}{4}}(1-{\frac {\hat {h}}{2}}))y_{n}$

$k_{4}={\hat {h}}(1-2{\hat {h}}+{\hat {h}}(1-{\frac {\hat {h}}{2}})+2{\hat {h}}(1+{\frac {3{\hat {h}}}{4}}-{\frac {\hat {h}}{4}}(1-{\frac {\hat {h}}{2}})))))y_{n}$

substitute these back in $y_{n+1}$, yields $y_{n+1}=[1+({\hat {h}})+{\frac {({\hat {h}})^{2}}{2}}+{\frac {({\hat {h}})^{3}}{6}}+{\frac {({\hat {h}})^{4}}{24}}]y_{n}=R({\hat {h}})y_{n}$

and so the characteristic polynomial

$P(z)=z-R(z);z={\hat {h}}$ for the absolute stability region for this method, set |R(z)|<1, and so we get the region in figure 1. case$2_{a}$

the table below shows the final forms for the stability function for different forms of RK4, these RK4's are different in the values of $b_{j}$, and they are fullfilling the consistencty requirement for the method i.e :$\sum _{j=1}^{i-1}a_{ij}=c_{i}\ \mathrm {for} \ i=2,\ldots ,s.$

case# $(b_{1},b_{2},b_{3},b_{4})$	stability function
caseI $({\frac {1}{8}},{\frac {3}{8}},{\frac {3}{8}},{\frac {1}{8}})$	$1+z+{\frac {z^{2}}{2}}+{\frac {z^{3}}{6}}+{\frac {z^{4}}{24}}$
caseII_a $({\frac {1}{6}},0,{\frac {2}{3}},{\frac {1}{6}})$	$1+z+{\frac {z^{2}}{2}}+{\frac {z^{3}}{6}}+{\frac {z^{4}}{24}}$
caseII_b $({\frac {1}{6}},0,{\frac {2}{3}},{\frac {1}{6}})$	$1+z+{\frac {5z^{2}}{6}}-{\frac {17z^{3}}{12}}-{\frac {2z^{4}}{24}}+{\frac {z^{5}}{12}}$
caseIII $({\frac {1}{12}},{\frac {2}{3}},{\frac {1}{12}},{\frac {1}{6}})$	$1+{\frac {11z}{12}}+{\frac {5z^{2}}{12}}+{\frac {z^{3}}{6}}+{\frac {z^{4}}{24}}$
caseIV $({\frac {1}{6}},0,{\frac {2}{3}},{\frac {1}{6}})$	$1+z+{\frac {z^{2}}{2}}+{\frac {z^{3}}{6}}+{\frac {z^{4}}{24}}$
classical Rk $({\frac {1}{6}},{\frac {1}{3}},{\frac {1}{3}},{\frac {1}{6}})$	$1+z+{\frac {z^{2}}{2}}+{\frac {z^{3}}{6}}+{\frac {z^{4}}{24}}$

see [2] to find the table

Plotting the stability region

In order to plot the stability region, we can set the stability function to be bounded by 1 and solve for the values of z, then draw z in the complex plane. Since R(z) is the unit circle in the complex plane, each point on the boundary can be represented as $e^{i\theta }$ and so by changing $\theta$ over the interval$[0,2\pi ]$, we can draw the boundaries of that region. The following w:OCTAVE/Matlab code does this by plotting contour curves until reaching the boudaries:

[x,y] = meshgrid(-6:0.01:6,-6:0.01:6);
z = x+i*y;
R = 1+z+0.5*z.^2+(1/6)*z.^3+(1/24)*z.^4;
zlevel = abs(R);
contour(x,y,zlevel,[1 1]);

The figure at right shows the absolute stability regions for RK4 cases which is tabulated above[3]

References

• Eberly, David (2008), stability analysis for systems of defferential equation.[4]
• Ababneh, Osama; Ahmad, Rokiah; Ismail, Eddie (2009), "on cases of fourth-order Runge-Kutta methods", European journal of scientific Research.[5]
• Mathews, John; Fink, Kurtis (1992), numerical methods using matlab.
• Stability of Runge-Kutta Methods [6]