Trying to get region of stability of RADUA IIA

44 Views Asked by At

The method can be turned from its Butcher table into this set of explicit expressions:

$$k_1 = f\bigg(t_n + \frac{1}{3}h, y_n + \frac{5}{12}hk_1 + \frac{-1}{12}hk_2\bigg)$$ $$k_2 = f\bigg(t_n + h, y_n + \frac{3}{4}hk_1 + \frac{1}{4}hk_2\bigg)$$

$$y_{n+1} = y_n + h\frac{3}{4}k_1 + h\frac{1}{4}k_2$$

This is not enough to find the region of stability, so we pick the test function $y' = \lambda y$ and apply the method to it. Shutting your brain down and doing mindless symbol substitution yields:

$$k_1 = \lambda(y_n + \frac{5}{12}hk_1 + \frac{-1}{12}hk_2)$$ $$k_2 = \lambda(y_n + \frac{3}{4}hk_1 + \frac{1}{4}hk_2)$$

$$k_1 = f\bigg(t_n + \frac{1}{3}h, y_n + \frac{5}{12}hk_1 + \frac{-1}{12}hk_2\bigg)$$ $$k_2 = f\bigg(t_n + h, y_n +\frac{3}{4}k_1 + \frac{1}{4}k_2\bigg)$$

$$y_{n+1} = y_n + h\frac{3}{4}k_1 + h\frac{1}{4}k_2$$

Applying this to the test function $y'=\lambda y$ yields:

$$k_1 = \lambda(y_n + \frac{5}{12}hk_1 + \frac{-1}{12}hk_2)$$ $$k_2 = \lambda(y_n + \frac{3}{4}hk_1 + \frac{1}{4}hk_2)$$

So:

$$k_1 = \frac{- 2 \lambda^{2} y_{n} + 6 \lambda y_{n}}{\lambda^{2} h - 4 \lambda h + 6 h}, \ k_2 = \frac{2 \lambda^{2} y_{n} + 6 \lambda y_{n}}{\lambda^{2} h - 4 \lambda h + 6 h}$$

This then yields:

$$y_{n+1} = y_n + h\frac{3}{4}\frac{- 2 \lambda^{2} + 6 \lambda}{\lambda^{2} h - 4 \lambda h + 6 h}y_{n} + h\frac{1}{4}\frac{2 \lambda^{2} + 6 \lambda}{\lambda^{2} h - 4 \lambda h + 6 h}y_{n}$$

So:

$$y_{n+1} = y_n\bigg( 1 + \frac{3}{4}\frac{- 2 \lambda^{2} + 6 \lambda}{\lambda^{2} - 4 \lambda + 6} + \frac{1}{4}\frac{2 \lambda^{2} + 6 \lambda}{\lambda^{2} - 4 \lambda + 6}\bigg)$$

Simplifying the above through sympy gives:

$$y_{n+1} = \frac{2 \left(\lambda + 3\right)}{\lambda^{2} - 4 \lambda + 6}y_{n}$$

We can rewrite the factor by substituting $\lambda$ with $z$ as:

$$R(z) = \frac{2z + 3}{z^2 - 4z + 6}$$

And we are now interested in $|R(z)| \leq 1$.

Solving the fully expanded resulting expression is proving time consuming and sympy is also struggling at finding an explicit expression. I am not sure if a mistake was made. What should I do to get an explicit representation of the region of stability?