Part 4: Stability Analysis - ASTR 596: Modeling the Universe

Learning Outcomes¶

By the end of this section, you will be able to:

Analyze stability regions for different integration methods
Diagnose and handle stiff equations
Predict when numerical methods will become unstable
Choose appropriate methods based on stability requirements
Implement implicit methods for challenging problems

Linear Stability Theory¶

To understand when methods fail catastrophically, we analyze their behavior on the test equation:

\frac{dy}{dt} = \lambda y

(1)

where $\lambda \in \mathbb{C}$ . The true solution is $y(t) = y_0 e^{\lambda t}$ .

When we apply a numerical method, we get:

y_{n+1} = R(h\lambda) y_n

(2)

where $R(z)$ is the stability function with $z = h\lambda$ .

For bounded solutions, we need $|R(h\lambda)| \leq 1$ . The stability region is the set of $z$ values where this holds.

Stability Functions for Common Methods¶

Let’s derive the stability functions for our methods:

Euler’s Method¶

Applying Euler to $y' = \lambda y$ :

y_{n+1} = y_n + h\lambda y_n = (1 + h\lambda)y_n

(3)

Therefore: $\boxed{R(z) = 1 + z}$

This is stable when $|1 + z| \leq 1$ , which defines a circle of radius 1 centered at $z = -1$ .

RK2 (Midpoint)¶

Following through the RK2 steps:

k_1 = \lambda y_n

(4)

k_2 = \lambda(y_n + \frac{h}{2}k_1) = \lambda y_n(1 + \frac{z}{2})

(5)

y_{n+1} = y_n + hk_2 = y_n(1 + z + \frac{z^2}{2})

(6)

Therefore: $\boxed{R(z) = 1 + z + \frac{z^2}{2}}$

RK4¶

Through similar analysis (but more algebra):

\boxed{R(z) = 1 + z + \frac{z^2}{2} + \frac{z^3}{6} + \frac{z^4}{24}}

(7)

This matches the Taylor series of $e^z$ through fourth order, but is not exactly $e^z$ truncated!

Physical Interpretation¶

For different types of problems, we need different stability characteristics:

Oscillatory Systems (Circular Orbits)¶

For a harmonic oscillator with frequency $\omega$ :

Eigenvalues: $\lambda = \pm i\omega$ (purely imaginary)
Stability requires the imaginary axis to be in the stability region
Euler: Marginally stable for $h < 2/\omega$ , but introduces artificial damping
RK4: Stable for $h < 2.8/\omega$
Leapfrog: Stable for $h < 2/\omega$ with no damping!

Dissipative Systems (Damped Motion)¶

For exponential decay with rate $\gamma$ :

Eigenvalue: $\lambda = -\gamma$ (negative real)
Stability requires the negative real axis in the stability region
Euler: Stable for $h < 2/\gamma$
RK4: Stable for $h < 2.78/\gamma$

Stiff Equations - When Stability Dominates¶

What Makes an Equation Stiff?¶

A stiff equation contains widely separated timescales. After fast transients die out, we’re left tracking slow evolution, but explicit methods still need tiny timesteps for stability.

Understanding Stiffness Through Examples¶

Example 1: Chemical Kinetics¶

Consider a reaction where species A slowly converts to B, but B rapidly equilibrates:

\frac{dA}{dt} = -0.04A

(8)

\frac{dB}{dt} = 0.04A - 10^4 B

(9)

The eigenvalues are $\lambda_1 = -0.04$ (slow) and $\lambda_2 = -10^4$ (fast).

After the initial transient ( $t > 0.001$ ), B has equilibrated and we’re just tracking the slow decay of A. But explicit methods need $h < 2/10^4 = 0.0002$ for stability, even though nothing interesting happens on this timescale!

Example 2: Stellar Interior with Radiative Cooling¶

In stellar shocks, radiative cooling can create extreme stiffness:

\frac{dT}{dt} = \text{heating} - \text{cooling}(T)

(10)

where cooling $\propto T^4$ for high temperatures. Near equilibrium:

Heating timescale: years
Cooling timescale: seconds
Stiffness ratio: 10⁷!

Non-Stiff vs Stiff: The Key Distinction¶

Non-stiff equations have eigenvalues with similar magnitudes. The timestep is chosen for accuracy, and stability is automatically satisfied. Examples:

Planetary orbits (all periods within a few orders of magnitude)
Pendulum oscillations
Wave propagation in uniform media

Stiff equations have eigenvalues with vastly different magnitudes. The timestep is forced by stability of the fastest mode, even after it has decayed. Examples:

Chemical reaction networks
Radiative transfer with cooling
Electrical circuits with multiple timescales

Implicit Methods for Stiff Problems¶

Implicit methods evaluate the derivative at the future point, giving superior stability at the cost of solving equations.

Backward Euler¶

Instead of $y_{n+1} = y_n + hf(y_n, t_n)$ , we use:

\boxed{y_{n+1} = y_n + hf(y_{n+1}, t_{n+1})}

(11)

This is implicit because $y_{n+1}$ appears on both sides!

Stability Analysis¶

For the test equation $y' = \lambda y$ :

y_{n+1} = y_n + h\lambda y_{n+1}

(12)

y_{n+1}(1 - h\lambda) = y_n

(13)

y_{n+1} = \frac{y_n}{1 - h\lambda}

(14)

Stability function: $\boxed{R(z) = \frac{1}{1-z}}$

This is stable for the entire left half-plane! We call this A-stable.

Solving the Implicit Equation¶

The challenge: we must solve $y_{n+1} = y_n + hf(y_{n+1}, t_{n+1})$ for $y_{n+1}$ .

For nonlinear $f$ , this requires iteration. Newton’s method is typical:

def backward_euler_step(y, t, h, f, df_dy, tol=1e-10):
    """
    Single backward Euler step
    Requires Jacobian df/dy for Newton iteration
    """
    y_new = y  # Initial guess
    
    for _ in range(10):  # Newton iterations
        # Residual: R(y_new) = y_new - y - h*f(y_new, t+h)
        residual = y_new - y - h*f(y_new, t + h)
        
        if np.linalg.norm(residual) < tol:
            return y_new
        
        # Jacobian of residual
        J = np.eye(len(y)) - h*df_dy(y_new, t + h)
        
        # Newton update
        delta = np.linalg.solve(J, -residual)
        y_new = y_new + delta
    
    raise RuntimeError("Newton iteration failed to converge")

Implicit Midpoint Rule¶

A second-order implicit method:

y_{n+1} = y_n + hf\left(\frac{y_n + y_{n+1}}{2}, t_n + \frac{h}{2}\right)

(15)

This is A-stable AND symplectic! It preserves geometric structure while handling stiffness.

The CFL Condition¶

The CFL condition is a necessary stability criterion for hyperbolic PDEs that also applies to ODEs from spatial discretization:

h \leq \frac{\Delta x}{c_{max}}

(16)

where $\Delta x$ is the spatial grid spacing and $c_{max}$ is the maximum wave speed.

For the wave equation $\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}$ discretized in space:

Eigenvalues: $\lambda \approx \pm ic/\Delta x$
Stability requires: $h < \Delta x/c$

This means finer spatial grids require smaller timesteps!

Diagnosing Integration Problems¶

Systematic Debugging Approach¶

When your integration fails, follow this diagnostic sequence:

Check for NaN/Inf: Immediate indicator of instability
Plot the solution: Look for oscillations growing exponentially
Compute eigenvalues: Estimate $\lambda_{max}$
Check stability limit: Is $h|\lambda_{max}| < \text{stability limit}$ ?
Test with smaller h: Does halving h fix the problem?
Monitor energy: For conservative systems, is energy bounded?

Stability Test Code¶

def stability_test(method, lambda_val, h_values):
    """Test stability for different timesteps"""
    results = []
    
    for h in h_values:
        y = 1.0  # Initial condition
        stable = True
        
        for n in range(1000):
            y = method(y, lambda_val, h)
            
            if abs(y) > 1e10:  # Explosion detected
                stable = False
                break
        
        results.append((h, stable, abs(y)))
    
    return results

Bridge to Part 5: From Stability to Speed¶

Understanding stability tells us what timesteps are allowed, but performance determines whether our simulations finish in reasonable time. Even the most stable method is useless if it takes years to run.

In Part 5, we’ll explore performance optimization through vectorization. Modern processors can perform multiple operations simultaneously, but only if we structure our code correctly. The difference between naive loops and optimized array operations can be 100× in execution speed—the difference between waiting hours and waiting weeks for results.

Next: Part 5 - Performance Optimization