Runge Kutta stability

2.1k Views Asked by At

I am facing a problem solving a ODE with a Runge-Kutta 4th order method:

The expression in order to solve is :

\begin{equation} Ay^{''}+By^{'}+Cy= Cu \end{equation}

\begin{equation} y =OUTPUT \end{equation}

\begin{equation} u=INPUT \end{equation}

We are using a sample time of 0.01s (with 0.001s sample time it does work), but the solver with some combinations of A,B and C does not converge. Then, we would like to know these combinations of A,B and C that make the method fail before start to using it.

Thanks in advance.

Case of study 1 (This crashes) $A=1.039\cdot10^3$ $B=5\cdot10^5$ $C=1.55\cdot10^7$

Case of study 2 (This works) $A=7.93\cdot10^3$ $B=5\cdot10^5$ $C=1.55\cdot10^5$

1

There are 1 best solutions below

1
On BEST ANSWER

Your problem is stiff. I recommend you to take a look at the book Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems by E. Hairer and G. Wanner. For stiff problems you need either an explicit method with small timestep or some stiff-stable implicit method.

To determine whether and RK method with timestep $\Delta t$ is applicable to some problem, you need to know two things:

  • Eigenvalues of your problem
  • Region of stability of the method used

The eigenvalues of your problems can be easily found from characteristic equation $$ A \lambda^2 + B\lambda + C = 0 $$ and that are $$ \lambda_1 \approx -447.927, \lambda_2 \approx -33.305 $$ for the unstable case and $$ \lambda_{1,2} \approx -31.5259 \pm 30.9955i $$ for the stable one.

For a method to be stable, all the values $z_i = \lambda_i \Delta t$ must lie in the method's stability region, which for RK4 is given by $$ \left|1 + z + \frac{z^2}{2} + \frac{z^3}{6} + \frac{z^4}{24}\right| \leq 1 $$ Here's an image of the stability region for RK4

RK4 stability region

You can see, that for $\Delta t = 0.001$ every $z_i$ lies in the stability region and for $\Delta t = 0.01$ the $\lambda_1 \approx -448$ gives $z_1 \approx -4.4$ which lies outside.