Skip to content

Fixed-Step Integrators

Fixed-step integrators use a constant time step throughout the integration. Unlike adaptive methods, they don't automatically adjust step size based on error estimates. They are simpler to implement and have predictable computational costs, but require careful step size selection to ensure accuracy. They provide regular output at fixed intervals, making them suitable for applications needing uniform sampling.

RK4: Classical Runge-Kutta

Brahe implements the classical 4th-order Runge-Kutta (RK4) method as it's primary fixed-step integrator. The 4th-order Runge-Kutta method (RK4) is the most popular fixed-step integrator, offering an excellent balance of accuracy and simplicity.

Algorithm

For \(\dot{\mathbf{x}} = \mathbf{f}(t, \mathbf{x})\), the RK4 method computes:

\[\begin{align} \mathbf{k}_1 &= \mathbf{f}(t, \mathbf{x}) \\ \mathbf{k}_2 &= \mathbf{f}(t + h/2, \mathbf{x} + h\mathbf{k}_1/2) \\ \mathbf{k}_3 &= \mathbf{f}(t + h/2, \mathbf{x} + h\mathbf{k}_2/2) \\ \mathbf{k}_4 &= \mathbf{f}(t + h, \mathbf{x} + h\mathbf{k}_3) \end{align}\]

The next state is then given by:

\[\mathbf{x}(t + h) = \mathbf{x}(t) + \frac{h}{6}(\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4)\]

Choosing Step Size

The step size h must balance accuracy and computational cost. Too large causes unacceptable errors; too small wastes computation. A decent starting point is to relate h to the characteristic time scale of the dynamics. For orbital dynamics, a common heuristic is

\[ h \approx \frac{T}{100 \text{ \textemdash } 1000} \]

where T is the orbital period.

Since fixed-step methods lack automatic error control, it is critical to validate that the step-size choice achieves the desired level of accuracy. Common validation approaches include:

  1. Analytical solution: Compare against closed-form solution (when available)
  2. Step Size Comparison: Run with both \(h\) and \(h/2\), compare results to confirm convergence
  3. Energy/momentum conservation: If you have a conserative system (in astrodynamics, this would be a gravitional-only system), check that total energy and angular momentum remain constant over time.
  4. Reference integrator: Compare against adaptive integrator with tight tolerances

Basic Integration Example

The following example demonstrates using the RK4 fixed-step integrator to integrate a simple harmonic oscillator.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
import brahe as bh
import numpy as np

# Define simple harmonic oscillator
omega = 1.0


def dynamics(t, state):
    x, v = state
    return np.array([v, -(omega**2) * x])


# Analytical solution
def analytical(t, x0=1.0, v0=0.0):
    x = x0 * np.cos(omega * t) + (v0 / omega) * np.sin(omega * t)
    v = -omega * x0 * np.sin(omega * t) + v0 * np.cos(omega * t)
    return np.array([x, v])


# Initial conditions
state0 = np.array([1.0, 0.0])
t_end = 4 * np.pi  # Two periods

print("RK4 Fixed-Step Integration Demonstration")
print("System: Simple Harmonic Oscillator (ω=1.0)")
print(f"Integration time: 0 to {t_end:.2f} (2 periods)\n")

# Test different step sizes
step_sizes = [0.5, 0.2, 0.1, 0.05]

for dt in step_sizes:
    config = bh.IntegratorConfig.fixed_step(step_size=dt)
    integrator = bh.RK4Integrator(2, dynamics, config=config)

    t = 0.0
    state = state0.copy()
    steps = 0

    # Integrate to end
    while t < t_end - 1e-10:
        state = integrator.step(t, state, dt)
        t += dt
        steps += 1

    # Compare with analytical solution
    exact = analytical(t)
    error = np.linalg.norm(state - exact)

    print(f"Step size dt={dt:5.2f}:")
    print(f"  Steps:      {steps}")
    print(f"  Final state: [{state[0]:.6f}, {state[1]:.6f}]")
    print(f"  Exact:       [{exact[0]:.6f}, {exact[1]:.6f}]")
    print(f"  Error:       {error:.2e}")
    print()

# Expected Output:
# RK4 Fixed-Step Integration Demonstration
# System: Simple Harmonic Oscillator (ω=1.0)
# Integration time: 0 to 12.57 (2 periods)

# Step size dt= 0.50:
#   Steps:      26
#   Final state: [0.907541, -0.413422]
#   Exact:       [0.907447, -0.420167]
#   Error:       6.75e-03

# Step size dt= 0.20:
#   Steps:      63
#   Final state: [0.999412, -0.033457]
#   Exact:       [0.999435, -0.033623]
#   Error:       1.68e-04

# Step size dt= 0.10:
#   Steps:      126
#   Final state: [0.999434, -0.033613]
#   Exact:       [0.999435, -0.033623]
#   Error:       1.05e-05

# Step size dt= 0.05:
#   Steps:      252
#   Final state: [0.999435, -0.033622]
#   Exact:       [0.999435, -0.033623]
#   Error:       6.56e-07

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
use brahe::integrators::*;
use nalgebra::DVector;

fn main() {
    // Define simple harmonic oscillator
    let omega: f64 = 1.0;
    let dynamics = move |_t: f64, state: &DVector<f64>, _params: Option<&DVector<f64>>| -> DVector<f64> {
        let x = state[0];
        let v = state[1];
        DVector::from_vec(vec![v, -omega.powi(2) * x])
    };

    // Analytical solution
    let analytical = |t: f64, x0: f64, v0: f64| -> (f64, f64) {
        let x = x0 * (omega * t).cos() + (v0 / omega) * (omega * t).sin();
        let v = -omega * x0 * (omega * t).sin() + v0 * (omega * t).cos();
        (x, v)
    };

    // Initial conditions
    let state0 = DVector::from_vec(vec![1.0, 0.0]);
    let t_end = 4.0 * std::f64::consts::PI;  // Two periods

    println!("RK4 Fixed-Step Integration Demonstration");
    println!("System: Simple Harmonic Oscillator (ω=1.0)");
    println!("Integration time: 0 to {:.2} (2 periods)\n", t_end);

    // Test different step sizes
    let step_sizes = vec![0.5, 0.2, 0.1, 0.05];

    for &dt in &step_sizes {
        let config = IntegratorConfig::fixed_step(dt);
        let integrator = RK4DIntegrator::with_config(2, Box::new(dynamics), None, None, None, config);

        let mut t = 0.0;
        let mut state = state0.clone();
        let mut steps = 0;

        // Integrate to end
        while t < t_end - 1e-10 {
            state = integrator.step(t, state, None, None).state;
            t += dt;
            steps += 1;
        }

        // Compare with analytical solution
        let (exact_x, exact_v) = analytical(t, 1.0, 0.0);
        let error = ((state[0] - exact_x).powi(2) + (state[1] - exact_v).powi(2)).sqrt();

        println!("Step size dt={:5.2}:", dt);
        println!("  Steps:      {}", steps);
        println!("  Final state: [{:.6}, {:.6}]", state[0], state[1]);
        println!("  Exact:       [{:.6}, {:.6}]", exact_x, exact_v);
        println!("  Error:       {:.2e}", error);
        println!();
    }
}

// Expected output:
// RK4 Fixed-Step Integration Demonstration
// System: Simple Harmonic Oscillator (ω=1.0)
// Integration time: 0 to 12.57 (2 periods)

// Step size dt= 0.50:
//   Steps:      26
//   Final state: [0.907541, -0.413422]
//   Exact:       [0.907447, -0.420167]
//   Error:       6.75e-03

// Step size dt= 0.20:
//   Steps:      63
//   Final state: [0.999412, -0.033457]
//   Exact:       [0.999435, -0.033623]
//   Error:       1.68e-04

// Step size dt= 0.10:
//   Steps:      126
//   Final state: [0.999434, -0.033613]
//   Exact:       [0.999435, -0.033623]
//   Error:       1.05e-05

// Step size dt= 0.05:
//   Steps:      252
//   Final state: [0.999435, -0.033622]
//   Exact:       [0.999435, -0.033623]
//   Error:       6.56e-07

See Also