Homework 1: ODE’s and Particle Dynamics

Homework #1: ODE’s and Particle Dynamics

Huai-Ping Lee

September 20, 2007

Implementation

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.)

Physical and Numerical Stability

Problem A. For a particle system that is only affected by gravity and drag force, we have the equations

[ ] [ ] [ ] ˙x = J x + vo , ¨x ˙x g

where J is the Jacobian matrix. The physical stability of the linear system depends on the eigenvalues of J: every solution is stable if all the eigenvalues have negative real part, or all the eigenvalues are less than or equal to zero if they are all real numbers [1]. In our case,

[ ] 0 1 J = 0 - kd∕m ,

and the eigenvalues are λ = 0, -k_d∕m. Therefore this linear system is physically stable as long as k_d ≥ 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

[ ] [ ] ˙x = J x ¨x ˙x ¨x = - ksx. m

We can then analyze the numerical stability of different methods using this equation. The first method we used is Euler step:

Xi+1 = Xi + h ˙Xi = Xi + h(JXi ) = (I + hJ)Xi,

where X =

. The result is an exponential function that converges only if the magnitude of eigenvalue of (I + hJ) is less than or equal to 1, i.e., |1 + hλ|≤ 1, where λ is the eigenvalue of J.

Another numerical method implemented is fourth-order Runge-Kutta. In this case,

k1 = hλXi k2 = hλ(1+ hλ ∕2)Xi k3 = hλ(1+ (hλ∕2)(1+ hλ∕2))Xi k4 = hλ(1+ hλ (1 + (hλ ∕2)(1+ hλ∕2))Xi ( 1 2 1 3 1 4) Xi+1 = 1 + hλ+ 2(h λ) + 6(hλ) + 24(hλ) Xi = q(hλ)Xi.

(Here we write λ in place of J for simplicity.) So we need |q(hλ)|≤ 1 to assure stability.

In our case,

[ 0 1 ] J = - ks- 0 , m

and λ = ± ∘ - ks∕m = ±i ∘ks-∕m- , where i = √ - 1 . Therefore |1 + hλ| = ∘1-+-h2(ks∕m-)- 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 k_s = 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 + kd-
m + ks
m x = 0, and the Jacobian matrix becomes

[ ] 0 1 J = - kms - kmd ,

so the eigenvalues are the roots of λ² + kmd λ + kms = 0, namely

∘ --------- - kmd± kmd2 - 4ksm λ = --------2-------.

Then there is a chance that |1 + hλ|≤ 1 (λ can be a complex number) so that the stability of explicit Euler method is assured.

Results

Here are some plots from our experiments. This is a case when Euler method blows up while RK4 is still stable.

Figure 1:

Euler method with k_s = 30, m = 3, and h = 0.1

Figure 2:

Fourth-order Runge-Kutta method with k_s = 30, m = 3, and h = 0.1

Summary

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.

References

[1] David Eberly. Stability analysis for systems of differential equations, February 2003.