Events and Multistep Solvers

Repeated discontinuities can change which stiff solver is most efficient. A multistep method such as BDF uses solution history to achieve high order. A state jump invalidates that history, so BDF must restart at first order. A one-step method such as ESDIRK34 has no multistep history to discard and can continue with its fixed third-order formula after an event. This may mean that BDF performs poorly in comparison with a one-step method, but this depends on many factors, even with the restart BDF might perform significantly better between events, so it is always worth trying both methods to see which is better.

The performance-solver-events example constructs a simple example that highlights the slowdown effect that can affect the BDF solver. It compares BDF with ESDIRK34 using dense nalgebra vectors and matrices and an LU linear solver:

type M = NalgebraMat<f64>;
type LS = NalgebraLU<f64>;

A Stiff Slow Manifold

The model is a stiff system with a slowly varying solution. We will perturb the system at regular events in such a way as to introduce discontinuities while maintaining a slowly varying solution (previous versions of this problem showed that the BDF solver did better at fast transients after the event, reducing the effect we are trying to demonstrate).

$$ \begin{aligned} y_1' &= -y_1, \\ y_2' &= 999y_1 - 1000y_2, \end{aligned} $$

with \(y(0)=(1,1)^T\). Its eigenvalues are \(-1\) and \(-1000\), giving a stiffness ratio of (1000). The vector ((1,1)^T) is the slow eigenvector.

At the event times

$$ t_k = 0.05k, \qquad k=1,\ldots,100, $$

the example applies the jump

$$ y(t_k^+) = y(t_k^-) + \begin{pmatrix}1\\1\end{pmatrix}. $$

Each jump stays on the slow eigenspace so we don't have any large fast transients. This is done so the benchmark is dominated by the cost of BDF's history restart, rather than the cost of repeatedly resolving the fast transients.

The model definition is:

    OdeBuilder::<M>::new()
        .rtol(RTOL)
        .atol([ATOL, ATOL])
        .rhs_implicit(
            |y, _p, _t, dy| {
                dy[0] = -LAMBDA_SLOW * y[0];
                dy[1] = (LAMBDA_FAST - LAMBDA_SLOW) * y[0] - LAMBDA_FAST * y[1];
            },
            |_y, _p, _t, v, jv| {
                jv[0] = -LAMBDA_SLOW * v[0];
                jv[1] = (LAMBDA_FAST - LAMBDA_SLOW) * v[0] - LAMBDA_FAST * v[1];
            },
        )
        .init(
            |_p, _t, y| {
                y[0] = 1.0;
                y[1] = 1.0;
            },
            2,
        )
        .build()
        .unwrap()

Restarting

After every event, the example updates both the state and its derivative. This keeps the solver state consistent before BDF reconstructs its first-order history. The loop stops exactly at each event, applies the jump, and then continues:

    let mut solver = problem.esdirk34::<LS>().unwrap();
    let started = Instant::now();

    for event_index in 1..NUM_EVENTS {
        solver
            .set_stop_time(EVENT_INTERVAL * event_index as f64)
            .unwrap();
        while !matches!(solver.step().unwrap(), OdeSolverStopReason::TstopReached) {}

        let state = solver.state_mut();
        state.y[0] += 1.0;
        state.y[1] += 1.0;
        state.dy[0] = -LAMBDA_SLOW * state.y[0];
        state.dy[1] = (LAMBDA_FAST - LAMBDA_SLOW) * state.y[0] - LAMBDA_FAST * state.y[1];
    }

    let elapsed = started.elapsed();
    let statistics = solver.get_statistics();
    #[cfg(test)]
    let state = solver.state();
    RunResult {
        solver: "ESDIRK34",
        elapsed,
        steps: statistics.number_of_steps,
        rejected_steps: statistics.number_of_error_test_failures,
        rhs_evals: solver.problem().eqn.statistics().number_of_calls,
        #[cfg(test)]
        final_state: [state.y[0], state.y[1]],
    }

    let mut solver = problem.bdf::<LS>().unwrap();
    let started = Instant::now();

    for event_index in 1..NUM_EVENTS {
        solver
            .set_stop_time(EVENT_INTERVAL * event_index as f64)
            .unwrap();
        while !matches!(solver.step().unwrap(), OdeSolverStopReason::TstopReached) {}

        let state = solver.state_mut();
        state.y[0] += 1.0;
        state.y[1] += 1.0;
        state.dy[0] = -LAMBDA_SLOW * state.y[0];
        state.dy[1] = (LAMBDA_FAST - LAMBDA_SLOW) * state.y[0] - LAMBDA_FAST * state.y[1];
    }

    let elapsed = started.elapsed();
    let statistics = solver.get_statistics();
    #[cfg(test)]
    let state = solver.state();
    RunResult {
        solver: "BDF",
        elapsed,
        steps: statistics.number_of_steps,
        rejected_steps: statistics.number_of_error_test_failures,
        rhs_evals: solver.problem().eqn.statistics().number_of_calls,
        #[cfg(test)]
        final_state: [state.y[0], state.y[1]],
    }

Results

The benchmark warms up both solvers and performs eleven timed runs of each. The example writes its most recent comparison to this table:

Wall-clock times depend on the machine, but this problem should show ESDIRK34 requiring substantially fewer steps than BDF at tight tolerances. BDF must repeatedly climb from order one after each jump instead of preserving high-order history over the full solve. However, it is not a general rule that one-step methods always beat BDF for event-driven stiff systems: repeatedly exciting a large fast mode can favour BDF instead, as we found for another version of this same problem.