## 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

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,
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

We can then analyze the numerical stability of different methods using this equation.
The first method we used is Euler step: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,

(Here we write λ in place of J for simplicity.) So we need |q(hλ)|≤ 1 to assure
stability.
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 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 + + 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.

### Results

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

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