Skip to content

Continuous

Fixed-step ODE integration plus a small library of canonical example systems used throughout the demos and tests.

All continuous models use SI units: - Displacement: metres (m) - Velocity: m/s - Angles: radians (rad)

Quick example

import numpy as np
import simweave as sw

msd = sw.MassSpringDamper(mass=1.0, damping=0.4, stiffness=4.0)
res = sw.simulate(msd, t_span=(0.0, 12.0), dt=0.01, x0=np.array([1.0, 0.0]))

Available systems

  • MassSpringDamper
  • SimplePendulum
  • QuarterCarModel
  • HalfCarModel
  • RollCarModel
  • FullCarModel
  • SeriesRLC
  • ThermalRC
  • TwoMassThermal

Each implements the DynamicSystem / SupportsDynamics protocol so custom systems plug into simulate() the same way.

Hybrid: continuous + discrete

Wrap a system in ContinuousProcess to register it on a SimEnvironment. The integrator advances one tick per env.tick() call, so continuous physics share the clock with queues and agents.

PID Control

Please see demo 20_PID_thermal_system.py to see a PID controller applied to a simple ThermalRC system. A temperature regulation system responding to workload changes could be helpful in real-world cases like: - CPUs - batteries - industrial heaters - HVAC systems

Vehicle dynamics examples

Quarter car

A 2 DOF (Degrees of Freedom) model with 4 variables (position + velocity).

  • z_s → sprung mass (body)
  • z_u → unsprung mass (wheel)

d = sw.QuarterCarModel(250, 40, 15000, 1500, 200000)
r = sw.simulate(d, (0.0, 2.0), dt=0.001, inputs=lambda t: 0.01)
See full demo: demos/10_quarter_car.py

Half car (pitch)

A 4 DOF (Degrees of Freedom) model for simulating pitch movement between front and back.

  • z_s → vertical body motion
  • theta → pitch
  • z_uf → front wheel
  • z_ur → rear wheel

model = sw.HalfCarModel(
    1200, 2500, 60, 60,
    20000, 20000,
    1500, 1500,
    150000, 150000,
    1.2, 1.6
)

r = sw.simulate(model, (0.0, 2.0), dt=0.001,
                inputs=lambda t: (0.01, 0.01))
See full demo: demos/17_half_car_pitch.py

Roll model (left/right)

A 4 DOF (Degrees of Freedom) model for simulating roll movement between left and right.

  • z_s → vertical body motion
  • phi → roll
  • z_ul → left wheel
  • z_ur → right wheel

model = sw.RollCarModel(
    1200, 2200, 60, 60,
    20000, 20000,
    1500, 1500,
    150000, 150000,
    1.6
)

r = sw.simulate(model, (0.0, 2.0), dt=0.001,
                inputs=lambda t: (0.01, 0.0))
See full demo: demos/18_half_car_roll_model.py

Full Car Model

A 7 DOF (Degrees of Freedom) model is included, extending prior examples.

  • z_s → vertical body motion
  • theta → pitch
  • phi → roll
  • z_ufl → front-left
  • z_ufr → front-right
  • z_url → rear-left
  • z_urr → rear-right

model = sw.FullCarModel(
    sprung_mass=1200,
    pitch_inertia=1500,
    roll_inertia=800,
    unsprung_mass=40,
    k_s=20000,
    c_s=1500,
    k_t=200000,
    a=1.2,
    b=1.3,
    track_width=1.6,
)

result = simulate(
    model, (0.0, 3.0),
    dt=0.001, 
    inputs=lambda t: np.array([0.01, 0.01, 0.0, 0.0]) if t > 1 else np.zeros(4),
)

fig = sw.plot_vehicle_metrics(
    res,
    title="Full car response (bump input)",
    model=model,
)
See full demo: demos/19_full_car_dynamics.py

Suspension Controllers

From v0.6.0, SimWeave includes controllers to help optimise dynamics according to prefered performance criteria: - Comfort (Acceleration of Vertical Vibrations) - Grip (Dynamic Tyre Loading) - Hybrid

These are realised by Skyhook and Groundhook control. See: - simweave.continuous.control.suspension.SkyhookDamper - simweave.continuous.control.suspension.GroundhookDamper - simweave.continuous.control.suspension.HybridActiveDamper

The former imagines one end of a damper 'connected to the sky' (an inertial reference point). The latter imagines one end affixed to the ground. This ensure dissipative force applies with respect to body movement, or wheel movement only as opposed to the releative velocities.

Example controller usage: 

passive = QuarterCarModel(250, 40, 15000, 1500, 200000)
controlled = QuarterCarModel(250, 40, 15000, 1500, 200000, controller=SkyhookDamper(1500))

r_passive = simulate(passive, (0, 2), dt=0.001, inputs=lambda t: 0.01)
r_control = simulate(controlled, (0, 2), dt=0.001, inputs=lambda t: 0.01)

z_passive = r_passive.state[:, 0]
z_control = r_control.state[:, 0]

print(f"Standard deviation controlled system: {z_control.std()}")
print(f"Standard deviation passive system: {z_passive.std()}")

A difference between an idealised behaviour, and the passive suspension can be used as a control signal. This could be used by an actuator in place of a damper, capable of adding energy into the system as a 'Fully Active' system. Alternatively, it could be a variable damping coefficient e.g. magnetorheological (MR) dampers.

To ensure fidelity with a real-world suspension that is only dissipative in nature, a wrapper allows for controller force to only apply when Force * relative velocity is negative.

See: simweave.continuous.control.suspension.SemiActiveWrapper

# enforce dissipative constraint
if F * v_rel > 0:
    return 0.0
Example wrapper usage: 
controller = SemiActiveWrapper(SkyhookDamper(1500))

API

Continuous-time dynamics: state-space solver, plug-in systems, hybrid wrapper.

DynamicSystem

Base class for first-order state-space dynamic systems.

SupportsDynamics

Bases: Protocol

Protocol for plug-in dynamic systems.

ContinuousProcess 

ContinuousProcess(system: SupportsDynamics, method: str = 'rk4', n_substeps: int = 1, inputs: InputFunc | None = None, name: str | None = None)

Bases: Entity

Drive a DynamicSystem from a SimEnvironment's tick loop.

Each tick, integrates n_substeps sub-intervals of length dt / n_substeps. The previous state is appended to history_state with the corresponding time in history_time.

Source code in src/simweave/continuous/solver.py 
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
def __init__(
    self,
    system: SupportsDynamics,
    method: str = "rk4",
    n_substeps: int = 1,
    inputs: InputFunc | None = None,
    name: str | None = None,
) -> None:
    super().__init__(
        name=name or getattr(system, "name", system.__class__.__name__)
    )
    if method not in _STEPPERS:
        raise ValueError(f"Unknown method '{method}'")
    if n_substeps < 1:
        raise ValueError("n_substeps must be >= 1")
    self.system = system
    self.method = method
    self.n_substeps = n_substeps
    self.inputs = inputs
    self.state = _as_state_vector(system.initial_state())
    self._stepper = _STEPPERS[method]
    self.history_time: list[float] = []
    self.history_state: list[np.ndarray] = []

PIDController 

PIDController(Kp: float, Ki: float = 0.0, Kd: float = 0.0, setpoint: float = 0.0)

Generic PID controller.

u(t) = Kp * e + Ki * ∫e dt + Kd * de/dt

Source code in src/simweave/continuous/control/pid.py 
11
12
13
14
15
16
17
18
19
20
21
22
23
24
def __init__(
    self,
    Kp: float,
    Ki: float = 0.0,
    Kd: float = 0.0,
    setpoint: float = 0.0,
):
    self.Kp = float(Kp)
    self.Ki = float(Ki)
    self.Kd = float(Kd)
    self.setpoint = float(setpoint)

    self._integral = 0.0
    self._prev_error = None

SkyhookDamper 

SkyhookDamper(damping: float)

Skyhook damping controller.

Applies force proportional to body velocity.

Source code in src/simweave/continuous/control/suspension.py 
7
8
def __init__(self, damping: float):
    self.c = float(damping)

GroundhookDamper 

GroundhookDamper(damping: float)

Groundhook damping controller.

Applies force proportional to wheel velocity.

Source code in src/simweave/continuous/control/suspension.py 
19
20
def __init__(self, damping: float):
    self.c = float(damping)

HybridActiveDamper 

HybridActiveDamper(c_sky: float, c_ground: float, alpha: float = 0.5)

Hybrid Skyhook & Groundhook.

Source code in src/simweave/continuous/control/suspension.py 
28
29
30
31
def __init__(self, c_sky: float, c_ground: float, alpha: float = 0.5):
    self.sky = SkyhookDamper(c_sky)
    self.ground = GroundhookDamper(c_ground)
    self.alpha = float(alpha)

SemiActiveWrapper 

SemiActiveWrapper(controller)

Enforces passive (dissipative) behaviour: F * v_rel <= 0.

Source code in src/simweave/continuous/control/suspension.py 
42
43
def __init__(self, controller):
    self.controller = controller

MassSpringDamper 

MassSpringDamper(mass: float | Mass, damping: float | Damping, stiffness: float | SpringStiffness, x0: tuple[float, float] = (0.0, 0.0))

Bases: DynamicSystem

Single degree-of-freedom mass-spring-damper.

State x = [displacement, velocity] Equation: m x'' + c x' + k x = F(t)

Source code in src/simweave/continuous/systems/mass_spring_damper.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
def __init__(
    self,
    mass: float | Mass,
    damping: float | Damping,
    stiffness: float | SpringStiffness,
    x0: tuple[float, float] = (0.0, 0.0),
):
    if mass <= 0:
        raise ValueError("mass must be positive")
    self.mass = self._val(mass)
    self.damping = self._val(damping)
    self.stiffness = self._val(stiffness)
    self._x0 = np.asarray([self._val(v) for v in x0], dtype=float)

SimplePendulum 

SimplePendulum(length: float | Distance, mass: float | Mass = 1.0, gravity: float | Acceleration = 9.81, damping: float | Damping = 0.0, x0: tuple[float, float] = (0.0, 0.0))

Bases: DynamicSystem

Simple pendulum with optional viscous damping and external torque.

State x = [theta, theta_dot] Equation: theta'' + (c / (m l^2)) theta' + (g/l) sin(theta) = tau / (m l^2)

Source code in src/simweave/continuous/systems/pendulum.py 
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
def __init__(
    self,
    length: float | Distance,
    mass: float | Mass = 1.0,
    gravity: float | Acceleration = 9.81,
    damping: float | Damping = 0.0,
    x0: tuple[float, float] = (0.0, 0.0),
):
    if length <= 0 or mass <= 0:
        raise ValueError("length and mass must be positive")
    self.length = self._val(length)
    self.mass = self._val(mass)
    self.gravity = self._val(gravity)
    self.damping = self._val(damping)
    self._x0 = np.asarray([self._val(v) for v in x0], dtype=float)
    self._inertia = self.mass * self.length**2

QuarterCarModel 

QuarterCarModel(sprung_mass: float | Mass, unsprung_mass: float | Mass, suspension_stiffness: float | SpringStiffness, damping: float | Damping, tyre_stiffness: float | SpringStiffness, x0: tuple[float, float, float, float] = (0.0, 0.0, 0.0, 0.0), controller=None)

Bases: DynamicSystem

Two degree-of-freedom quarter-car suspension model.

State vector x = [z_s, z_s_dot, z_u, z_u_dot] where: z_s: sprung mass displacement z_u: unsprung mass displacement

Source code in src/simweave/continuous/systems/quarter_car.py 
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
def __init__(
    self,
    sprung_mass: float | Mass,
    unsprung_mass: float | Mass,
    suspension_stiffness: float | SpringStiffness,
    damping: float | Damping,
    tyre_stiffness: float | SpringStiffness,
    x0: tuple[float, float, float, float] = (0.0, 0.0, 0.0, 0.0),
    controller = None
):
    if self._val(sprung_mass) <= 0 or self._val(unsprung_mass) <= 0:
        raise ValueError("Masses must be positive")
    self.m_s = self._val(sprung_mass)
    self.m_u = self._val(unsprung_mass)
    self.k_s = self._val(suspension_stiffness)
    self.c_s = self._val(damping)
    self.k_t = self._val(tyre_stiffness)
    self._x0 = np.asarray([self._val(v) for v in x0], dtype=float)
    self.controller = controller

HalfCarModel 

HalfCarModel(sprung_mass: float | Mass, pitch_inertia: float | Inertia, unsprung_mass_front: float | Mass, unsprung_mass_rear: float | Mass, k_sf: float | SpringStiffness, k_sr: float | SpringStiffness, c_sf: float | Damping, c_sr: float | Damping, k_tf: float | SpringStiffness, k_tr: float | SpringStiffness, a: float | Distance, b: float | Distance, x0: tuple[float, ...] = (0.0,) * 8, controller=None)

Bases: DynamicSystem

Half-car suspension model with pitch dynamics.

State vector: x = [z_s, z_s_dot, theta, theta_dot, z_uf, z_uf_dot, z_ur, z_ur_dot]

Source code in src/simweave/continuous/systems/half_car.py 
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
def __init__(
    self,
    sprung_mass: float | Mass,
    pitch_inertia: float | Inertia,
    unsprung_mass_front: float | Mass,
    unsprung_mass_rear: float | Mass,
    k_sf: float | SpringStiffness,
    k_sr: float | SpringStiffness,
    c_sf: float | Damping,
    c_sr: float | Damping,
    k_tf: float | SpringStiffness,
    k_tr: float | SpringStiffness,
    a: float | Distance,  # CG → front axle
    b: float | Distance,  # CG → rear axle
    x0: tuple[float, ...] = (0.0,) * 8,
    controller = None
):
    if self._val(sprung_mass) <= 0 or self._val(pitch_inertia) <= 0:
        raise ValueError("Mass and inertia must be positive")

    self.m_s = self._val(sprung_mass)
    self.I_y = self._val(pitch_inertia)

    self.m_uf = self._val(unsprung_mass_front)
    self.m_ur = self._val(unsprung_mass_rear)

    self.k_sf = self._val(k_sf)
    self.k_sr = self._val(k_sr)

    self.c_sf = self._val(c_sf)
    self.c_sr = self._val(c_sr)

    self.k_tf = self._val(k_tf)
    self.k_tr = self._val(k_tr)

    self.a = self._val(a)
    self.b = self._val(b)

    self._x0 = np.asarray([self._val(v) for v in x0], dtype=float)
    self.controller = controller

RollCarModel 

RollCarModel(sprung_mass: float | Mass, roll_inertia: float | Inertia, unsprung_mass_left: float | Mass, unsprung_mass_right: float | Mass, k_sl: float | SpringStiffness, k_sr: float | SpringStiffness, c_sl: float | Damping, c_sr: float | Damping, k_tl: float | SpringStiffness, k_tr: float | SpringStiffness, track_width: float | Distance, x0: tuple[float, ...] = (0.0,) * 8, controller=None)

Bases: DynamicSystem

Roll-only vehicle model (left-right dynamics).

Source code in src/simweave/continuous/systems/roll_car.py 
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
def __init__(
    self,
    sprung_mass: float | Mass,
    roll_inertia: float | Inertia,
    unsprung_mass_left: float | Mass,
    unsprung_mass_right: float | Mass,
    k_sl: float | SpringStiffness,
    k_sr: float | SpringStiffness,
    c_sl: float | Damping,
    c_sr: float | Damping,
    k_tl: float | SpringStiffness,
    k_tr: float | SpringStiffness,
    track_width: float | Distance,
    x0: tuple[float, ...] = (0.0,) * 8,
    controller = None,
):
    self.m_s = self._val(sprung_mass)
    self.I_phi = self._val(roll_inertia)

    self.m_ul = self._val(unsprung_mass_left)
    self.m_ur = self._val(unsprung_mass_right)

    self.k_sl = self._val(k_sl)
    self.k_sr = self._val(k_sr)

    self.c_sl = self._val(c_sl)
    self.c_sr = self._val(c_sr)

    self.k_tl = self._val(k_tl)
    self.k_tr = self._val(k_tr)

    self.w = self._val(track_width)

    self._x0 = np.asarray([self._val(v) for v in x0], dtype=float)
    self.controller = controller

FullCarModel 

FullCarModel(sprung_mass: float | Mass, pitch_inertia: float | Inertia, roll_inertia: float | Inertia, unsprung_mass: float | Mass, k_s: float | SpringStiffness, c_s: float | Damping, k_t: float | SpringStiffness, a: float | Distance, b: float | Distance, track_width: float | Distance, x0: tuple[float, ...] = (0.0,) * 14, controller=None)

Bases: DynamicSystem

Full car model with pitch and roll dynamics.

Source code in src/simweave/continuous/systems/full_car.py 
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
def __init__(
    self,
    sprung_mass: float | Mass,
    pitch_inertia: float | Inertia,
    roll_inertia: float | Inertia,
    unsprung_mass: float | Mass,
    k_s: float | SpringStiffness,
    c_s: float | Damping,
    k_t: float | SpringStiffness,
    a: float | Distance,       # CG → front axle
    b: float | Distance,       # CG → rear axle
    track_width: float | Distance,
    x0: tuple[float, ...] = (0.0,) * 14,
    controller = None,
):
    self.m_s = self._val(sprung_mass)
    self.I_y = self._val(pitch_inertia)
    self.I_phi = self._val(roll_inertia)

    self.m_u = self._val(unsprung_mass)

    self.k_s = self._val(k_s)
    self.c_s = self._val(c_s)
    self.k_t = self._val(k_t)

    self.a = self._val(a)
    self.b = self._val(b)
    self.w = self._val(track_width)

    self._x0 = np.asarray([self._val(v) for v in x0], dtype=float)
    self.controller = controller

SeriesRLC 

SeriesRLC(resistance: float | Resistance, inductance: float | Inductance, capacitance: float | Capacitance, x0: tuple[float, float] = (0.0, 0.0))

Bases: DynamicSystem

Series RLC circuit in charge-current form.

State x = [q, i] Equation: L q'' + R q' + (1/C) q = V(t)

Source code in src/simweave/continuous/systems/rlc.py 
21
22
23
24
25
26
27
28
29
30
31
32
33
def __init__(
    self,
    resistance: float | Resistance,
    inductance: float | Inductance,
    capacitance: float | Capacitance,
    x0: tuple[float, float] = (0.0, 0.0),
):
    if self._val(inductance) <= 0 or self._val(capacitance) <= 0:
        raise ValueError("inductance and capacitance must be positive")
    self.resistance = self._val(resistance)
    self.inductance = self._val(inductance)
    self.capacitance = self._val(capacitance)
    self._x0 = np.asarray([self._val(v) for v in x0], dtype=float)

ThermalRC 

ThermalRC(thermal_resistance: float | ThermalResistance, thermal_capacitance: float | ThermalCapacitance, ambient_temperature: float | Temperature = 293.15, initial_temperature: float | Temperature = 293.15, controller: None | PIDController = None)

Bases: DynamicSystem

First-order lumped-capacitance thermal model.

Parameters:

Name Type Description Default
thermal_resistance float | ThermalResistance

R_th in K/W. Larger => slower equilibration with ambient.

 required
thermal_capacitance float | ThermalCapacitance

C_th in J/K. Larger => more thermal inertia.

 required
ambient_temperature float | Temperature

T_inf in Kelvin (or any consistent unit). Default 293.15 K (20 C).

 293.15
initial_temperature float | Temperature

Starting body temperature, same unit as ambient_temperature.

 293.15
Source code in src/simweave/continuous/systems/thermal.py 
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
def __init__(
    self,
    thermal_resistance: float | ThermalResistance,
    thermal_capacitance: float | ThermalCapacitance,
    ambient_temperature: float | Temperature = 293.15,
    initial_temperature: float | Temperature = 293.15,
    controller: None | PIDController = None,
) -> None:
    if self._val(thermal_resistance) <= 0 or self._val(thermal_capacitance) <= 0:
        raise ValueError(
            "thermal_resistance and thermal_capacitance must be positive."
        )
    self.R_th = self._val(thermal_resistance)
    self.C_th = self._val(thermal_capacitance)
    self.T_inf = self._val(ambient_temperature)
    self._x0 = np.array([self._val(initial_temperature)], dtype=float)
    self.controller = controller

time_constant property 

time_constant: float

First-order time constant tau = R_th * C_th (seconds).

TwoMassThermal 

TwoMassThermal(C_core: float | ThermalCapacitance, C_sink: float | ThermalCapacitance, k_core_to_sink: float | ThermalConductance, R_sink_to_ambient: float | ThermalResistance, ambient_temperature: float | Temperature = 293.15, initial_core: float | Temperature = 293.15, initial_sink: float | Temperature = 293.15)

Bases: DynamicSystem

Two-lump thermal model: core <-> sink <-> ambient.

State: [T_core, T_sink], both in the same unit as ambient. Q_in(t) is applied to the core. Sink dissipates to ambient via R_sink_to_ambient; core couples to sink via conductance k_core_to_sink (W/K).

Source code in src/simweave/continuous/systems/thermal.py 
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
def __init__(
    self,
    C_core: float | ThermalCapacitance,
    C_sink: float | ThermalCapacitance,
    k_core_to_sink: float | ThermalConductance,
    R_sink_to_ambient: float | ThermalResistance,
    ambient_temperature: float | Temperature = 293.15,
    initial_core: float | Temperature = 293.15,
    initial_sink: float | Temperature = 293.15,
) -> None:
    for name, val in (
        ("C_core", C_core),
        ("C_sink", C_sink),
        ("k_core_to_sink", k_core_to_sink),
        ("R_sink_to_ambient", R_sink_to_ambient),
    ):
        if self._val(val) <= 0:
            raise ValueError(f"{name} must be positive.")
    self.C_core = self._val(C_core)
    self.C_sink = self._val(C_sink)
    self.k_cs = self._val(k_core_to_sink)
    self.R_sa = self._val(R_sink_to_ambient)
    self.T_inf = self._val(ambient_temperature)
    self._x0 = np.array([self._val(initial_core), self._val(initial_sink)], dtype=float)

simulate 

simulate(system: SupportsDynamics, t_span: tuple[float, float], dt: float, x0: Iterable[float] | ndarray | None = None, method: str = 'rk4', inputs: Optional[InputFunc] = None) -> SimulationResult

Integrate a system over t_span with fixed step dt.

Source code in src/simweave/continuous/solver.py 
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
def simulate(
    system: SupportsDynamics,
    t_span: tuple[float, float],
    dt: float,
    x0: Iterable[float] | np.ndarray | None = None,
    method: str = "rk4",
    inputs: Optional[InputFunc] = None,
) -> SimulationResult:
    """Integrate a system over ``t_span`` with fixed step ``dt``."""
    t0, tf = t_span
    if dt <= 0:
        raise ValueError("dt must be positive.")
    if tf <= t0:
        raise ValueError("t_span must satisfy tf > t0.")
    if method not in _STEPPERS:
        raise ValueError(f"Unknown method '{method}'. Use one of {tuple(_STEPPERS)}")

    x_init = _as_state_vector(system.initial_state() if x0 is None else x0)
    if hasattr(system, "validate_state"):
        system.validate_state(x_init)  # type: ignore[attr-defined]

    n_steps = int(np.floor((tf - t0) / dt)) + 1
    time = t0 + np.arange(n_steps + 1, dtype=float) * dt
    time = time[time <= tf + 1e-12]

    state = np.zeros((time.size, x_init.size), dtype=float)
    state[0] = x_init

    t0_input = _resolve_inputs(system, inputs, time[0])
    inputs_log = [np.atleast_1d(t0_input).astype(float) if t0_input is not None else np.zeros(1)]

    stepper = _STEPPERS[method]
    for i in range(1, time.size):
        t_prev = time[i - 1]

        # --- resolve inputs ---
        if inputs is not None:
            u = inputs(t_prev)
        else:
            u = system.inputs(t_prev)

        # normalise to array for storage
        if u is None:
            u_arr = np.zeros(1)
        else:
            u_arr = np.atleast_1d(u).astype(float)

        inputs_log.append(u_arr)

        # --- step ---
        state[i] = stepper(system, time[i - 1], state[i - 1], dt, inputs)

    # gather inputs
    inputs_array = np.vstack(inputs_log)

    labels = (
        tuple(system.state_labels())
        if hasattr(system, "state_labels")
        else tuple(f"x{i}" for i in range(x_init.size))
    )
    return SimulationResult(
        time=time,
        state=state,
        state_labels=labels,
        system_name=getattr(system, "name", system.__class__.__name__),
        method=method,
        inputs=inputs_array,
    )