Which Simulation Steps Need More Work?

Botie generates sequences of cooling temperatures for a learned simulator. A short fixed interval gives useful detail during a rapid change, but repeats unnecessary calculations once the temperature barely moves. A long interval is cheaper and can miss the part the model needed to learn.

We need a way to decide whether a proposed step is accurate enough before accepting it. One inspectable approach evaluates the same interval in two different ways. If their disagreement is too large for our stated tolerance, try a shorter interval from the same starting state.

Botie’s fixed-step simulator needs to identify intervals that are too coarse before using their states as training targets.

{
  "candidate": 9.5062,
  "estimated_error": 0.0062,
  "tolerance": 0.1951,
  "accepted": 1
}

A hot pump tends to cool toward the surrounding air. Our assumed model makes the cooling rate larger when the temperature difference is larger. The state we follow is that difference, called excess temperature above ambient.

The nonnegative rate parameter k controls how quickly the difference decays. A timestep h says how much physical time a proposed update would cover. Euler’s method uses the rate at the start of that interval to estimate the next state.

For now, keep an attempt separate from an accepted advance. If the later check rejects the approximation, the simulated clock and starting state must stay where they were.

rate_of_change = -k × excess_temperature
Euler(y,h) = y + h × rate_of_change

We are changing how the cooling rule is evaluated. Measurements are still needed to establish whether that rule describes the equipment.

Calculate one full Euler step from the starting state. Separately calculate a half step, evaluate the rate again at that midpoint, and take another half step. The second half must begin at the midpoint result.

Two half steps are not two copies of the first half-step value. Recomputing the rate is what gives the refined route additional information about the changing state.

coarse = Euler(y,h)
midpoint = Euler(y,h/2)
fine = Euler(midpoint,h/2)

y = 10, k = 0.5, h = 1
coarse = 5
midpoint = 7.5
fine = 5.625

Comparing a full step with only the first half step would compare different physical times.

temperature = 10.0
cooling_rate = 0.5
half_step = 0.5
midpoint = temperature - half_step * cooling_rate * temperature
endpoint = midpoint - half_step * cooling_rate * midpoint
print(midpoint)
print(endpoint)

7.5
5.625

The two routes used different amounts of detail to reach the same future time. Their disagreement gives us information about the approximation. For smooth behaviour at a sufficiently small step, Euler’s error order lets us use the magnitude of fine-minus-coarse as an estimate of the refined step’s leading local error.

The estimate still has assumptions. At a large step or across a sudden change in the rate, agreement can be misleading. Later we should compare whole trajectories with a reference solution or with a finer independent calculation.

Euler refined local error estimate ≈ abs(fine - coarse)
this attempt: abs(5.625 - 5) = 0.625

The error-order factor belongs to the method. This Euler expression should not be copied unchanged into a checker for an arbitrary solver.

An absolute tolerance is an allowance in the state’s units. A relative tolerance adds an allowance proportional to the candidate’s magnitude. Combining them lets an application express a useful scale near zero as well as away from it.

Choose tolerances from the downstream question. A display rounding choice is not evidence that a learning target can tolerate the same error. Multiple state variables may need separate absolute tolerances because their units and useful resolutions differ.

allowed local error = atol + rtol × abs(fine)
atol = 0.1, rtol = 0.01
allowed = 0.1 + 0.01 × 5.625 = 0.15625
0.625 exceeds 0.15625, so reject this attempt

Our teaching rule compares scalar values. Production solvers can combine scaled component errors using a method-specific norm.

On rejection, keep the starting state and starting time. Shorten the proposed interval and recompute both routes. Accepting the already-rejected state and merely shortening the next step would preserve the error we meant to avoid.

The programming challenge halves the interval until it passes, with a bound on the number of rejections. This simple controller is easy to inspect. Production controllers usually use the method order, error estimate and bounded scale changes to choose steps more efficiently.

A rejection bound or an inability to make a representable advance must produce an explicit failure, not an infinite retry loop.

Accepted local errors can accumulate over a long trajectory. Instability and a wrong physical model are separate sources of failure. Compare complete trajectories against an analytic answer where available, or against independent and refined numerical calculations.

For the stated constant-rate cooling model, the exponential solution provides a useful reference. Agreement with that solution verifies the numerical treatment of this model, not whether a real pump follows the same law.

exact excess at time t = initial_excess × exp(-k × t)

The learned simulator’s training labels inherit the numerical procedure and the physical assumptions that produced them.

A system can contain rapidly decaying behaviour alongside a much slower change we care about. An explicit method may need very short steps to keep the fast part stable even after it looks small. This is one reason a differential equation can be called stiff.

More function calls alone do not prove stiffness, but persistent difficulty is a reason to inspect the model and solver choice. SciPy offers implicit methods such as Radau and BDF for stiff problems. Replacing the solver still requires checking the resulting trajectory and diagnostics.

The workshop demonstrates local comparison and rejection only. It is not an implementation of RK45 or a general-purpose stiff solver.

In a real solver call, inspect success and the diagnostic message before treating the returned values as training data. Record the method, tolerances and initial state with the generated trajectory.

Requested output timestamps are also different from internal solver steps. A solver can take many adaptive internal advances and interpolate its result at the requested times. A small output table does not imply that only a few evaluations were performed.

Keep numerical settings in dataset provenance so that a later simulator revision does not silently change what a trained model learned to reproduce.

from scipy.integrate import solve_ivp

result = solve_ivp(
    lambda t,y: -0.5*y,
    (0.0,1.0), [10.0],
    method="RK45", rtol=1e-7, atol=1e-9,
    t_eval=[0.0,0.5,1.0],
)
if not result.success:
    raise RuntimeError(result.message)
print(result.t)
print(result.y[0])
print("function evaluations",result.nfev)