For both problems, we implemented two kinds of integrators: Euler method and 4th order Runge-Kutta method. The code can be found at ../code/PBM.zip. (This is a newer version, fixed some bugs in the demo version.)
Problem A. For a particle system that is only affected by gravity and drag force, we have the equations
and the eigenvalues are λ = 0, -kd∕m. Therefore this linear system is physically stable as long as kd ≥ 0. In this case there is no equilibrium state, so we will not discuss its numerical stability.
Problem B. For spring-mass systems, the equations become
Another numerical method implemented is fourth-order Runge-Kutta. In this case,
In our case,
and λ = ± = ±i, where i = . Therefore |1 + hλ| = is always greater than 1, and explicit Euler method is always unstable if there is no damping (!!) For RK4 method, it is hard to derive an inequality, but we can verify our example situation. Suppose we have ks = 30, m = 3, and h = 0.1, then we can calculate |q(hλ)|≈ 0.999993142337597 < 1. So the RK4 integrator is still numerically stable under this setting.
In general, if the equation of have some damping term dependent on , we will then have the equation + + x = 0, and the Jacobian matrix becomes
so the eigenvalues are the roots of λ2 + λ + = 0, namely
Then there is a chance that |1 + hλ|≤ 1 (λ can be a complex number) so that the stability of explicit Euler method is assured.
Here are some plots from our experiments. This is a case when Euler method blows up while RK4 is still stable.
RK4 is superior to Euler method in terms of accuracy and stability. However, it is not free. Since the RK4 method evaluates the integral four times to compute the result, it takes at least four times as much work as an Euler integrator.