Solver mathematics — the forms, the discretizations and the stepping scheme#

This page writes out the equations ICoreBlocks actually computes: the state-space and transfer-function forms a block holds, how a transfer function becomes a state space, the seven continuous-to-discrete conversions and the exact matrices each one produces, the integration schemes the simulator steps a continuous block with, and what step / impulse / ramp in the command window compute. Every formula is the one in the cited source file; nothing here is a textbook restatement. Sample-time inheritance and the one-sample loop delay are described on Sample time and loops — how the simulator paces a diagram; tolerances, the checks against runaway values and the fine print of each method's failure modes are on Numerics — what the solver will and will not do.

The two forms#

State space (src/ICoreBlocks/ICoreMath/ControlSystems/ICoreStateSpace.cpp). Every linear block holds six matrices, not four. The input is split into u, the block's ordinary input signal, and f, a second input channel meant to carry a nonlinear term that is computed outside the block and wired in as a signal:

continuous:   dx/dt = A x(t) + Bu u(t) + Bf f(t)
              y(t)  = C x(t) + Du u(t) + Df f(t)

discrete:     x[k+1] = A x[k] + Bu u[k] + Bf f[k]
              y[k]   = C x[k] + Du u[k] + Df f[k]

with n states, m_u inputs, m_f nonlinearity channels and p outputs. The dimension check (verifyMatrixDimensions) requires A square, Bu and Bf to have n rows, C to have n columns, Du/Df to have p rows and as many columns as Bu/Bf. A state space built from four matrices (A, B, C, D) gets Bf and Df as zero columns of the right height, so the plain linear blocks — State Space — Control Systems/Continuous, Discrete State Space — Control Systems/Discrete — are the special case f = 0. The Nonlinear State Space — Control Systems/Continuous block is the one whose two input ports are u and f, and whose config keys are literally A, Bu, Bf, C, Du, Df; its file banner says the same thing — the block is still linear in (u, f), and f is whatever the diagram feeds it. Ts <= 0 marks a state space as continuous (isContinuous()); Ts > 0 is a sampling period in seconds. Time never appears explicitly in A..Df — a time-varying block rewrites its matrices itself.

Transfer function (ICoreTransferFunction.cpp, Foundation/ICorePolynomial.h). One numerator and one denominator polynomial, coefficients written highest power first — tf([1],[1 2]) is 1 / (s + 2), and polynomial([..]) in the console is documented as "descending coefficients" (Command glossary — console commands, verbs, functions). The order of the transfer function is the order of its denominator (getOrder()). A numerator of higher order than the denominator — an improper model such as tf([1 2 3],[1 2]) — is held as a value, as in MATLAB; what refuses it is each operation that needs a proper model: a state-space realization ("an improper transfer function has no explicit state-space realization") and a time response ("Cannot simulate the time response of improper (non-causal) models."). As with a state space, Ts <= 0 is continuous (variable s), Ts > 0 discrete (variable z).

How a transfer function is realized as a state space#

ICoreStateSpace::fromTransferFunction() builds the controllable canonical form, the same realization MATLAB's tf2ss produces. With the denominator normalized so its leading coefficient is 1, den = s^n + a1 s^(n-1) + … + an, and the numerator zero-padded to the same length, num = b0 s^n + b1 s^(n-1) + … + bn:

A = [ -a1  -a2  …  -a(n-1)  -an ]      B = [ 1 ]
    [  1    0   …    0        0  ]          [ 0 ]
    [  0    1   …    0        0  ]          [ … ]
    [  …              1        0  ]          [ 0 ]

C = [ b1 - b0 a1,  b2 - b0 a2,  …,  bn - b0 an ]      D = [ b0 ]

The realized state space carries the transfer function's Ts. An order of 0 ("Invalid transfer function"), a zero (or < 1e-15) leading denominator coefficient ("Invalid transfer function denominator") or an improper model is refused, and an empty state space comes back. tf2ss(G) in the command window is this function; the run below shows tf([1 3],[1 2 5]) become A = [-2 -5; 1 0], B = [1; 0], C = [1 3], D = [0].

The seven discretization methods#

ICoreStateSpaceDiscretization::discretizeStateSpace(ss, Ts, method) (src/ICoreBlocks/ICoreMath/ControlSystems/StaticMethods/ICoreStateSpaceDiscretization.cpp) is the one place a continuous state space becomes a discrete one, whatever the caller: the simulator running in Discrete mode, the code exporters (targets always run the discrete model) and step/impulse/ramp on a continuous transfer function all go through it. Before the method runs, Bu|Bf are joined into one B and Du|Df into one D; afterwards the results are split back at column m_u. The continuous input must have Ts <= 0 and the requested Ts must be > 0, or the call logs an error and returns an empty state space. I is the n×n identity, Ts the sample period, and expm is the matrix exponential (ICoreMatrix::expm() in Foundation/ICoreMatrix.cpp, which delegates to Eigen's MatrixBase::exp()).

The strings in the table are what the discretizer matches, case-insensitively, after lower-casing ("first-order hold", "first order hold", "firstorderhold", "backward euler", "backwardeuler" and the forward variants are accepted too). Users do not type them: the solver setting Global Discretization Method (Solver Configuration panel, and getModelConfig discretizationMethod / setModelConfig discretizationMethod <name> in the console) takes one of the seven product names Zero-order Hold, First-order Hold, Impulse, Tustin, Matched, Backward Euler, Forward Euler (ICoreModelConfigurator.cpp), and each block's environment maps that to the short string when the model is built in Discrete mode (ICoreBlockSolverEnvironment::discretize(), src/ICoreBlocks/ICoreModel/SolverEnvironments/). The Ts handed to the discretizer is the block's own resolved sample time — see Sample time and loops — how the simulator paces a diagram.

An unrecognised method becomes Zero-order Hold, silently. Both dispatch points fall through to ZOH on no match — the else in discretize() and the else in discretizeStateSpace() — and neither logs. setModelConfig in the console refuses a name that is not one of the seven ("no option matches 'Tustn' — available: …", measured 2026-08-17), and the panel offers only the list, so the reachable path is a hand-edited project solver.ini ([Solver] discretizationMethod=…), which is read without validation despite the comment beside it (ICoreStudioSerialization.cpp; ICoreModelConfigurator::setDiscretizationMethod assigns the string as-is). A model discretized that way is valid, plausible and wrong. What each method does at the edges — a singular A under Tustin or Backward Euler, a Ts that aliases a mode — is on Numerics — what the solver will and will not do.

The stepping scheme#

The simulator (src/ICoreBlocks/ICoreSimulation/Core/ICoreModelSimulator.cpp) is a fixed-order, time-marching solver: at each solver time tn it visits every block in the build's solve order and asks it to bring its state and outputs to tn (ICoreBlockSolverEnvironment::solve(tn)). Two solver settings decide what that means:

• Solver — Continuous or Discrete. Under Discrete, every block that holds a continuous state space has it discretized once at build time (method above, at the block's Ts) and then runs the recursion x[k+1] = Ad x[k] + Bd u[k], y[k] = Cd x[k] + Dd u[k]: the block keeps x[k+1] from the previous visit, computes the output from the current state and input, then stores the next state. Under Continuous, blocks with continuous dynamics are integrated numerically (below); blocks that are discrete by nature (unit delays, discrete-time integrators, the discrete filters) run their recursion at their own rate in both modes.
• Stepping — which integrator, and whether the step is fixed. The choices, in the exact strings the setting holds (getModelConfig steppingType): Fixed-step | Euler | RK1, Fixed-step | Improved Euler (Heun) | RK2, Fixed-step | Kutta 3rd-order | RK3, Fixed-step | Runge-Kutta | RK4 (the default), Fixed-step | Backward Euler (implicit) | BE1, Fixed-step | Trapezoidal (implicit) | TR2, Variable-step | Dormand-Prince | RK45, Variable-step | Bogacki–Shampine | RK23, Variable-step | TR-BDF2 (implicit) | TRBDF2. A Discrete solver is always fixed-step. (RK3 was named Fixed-step | Bogacki–Shampine | RK3 until 2026-09-02; a solver.ini carrying the old string still selects it — setSteppingType.)

The first visit is the initial condition. At startTime every block is visited once and reports its state as it stands — an Integrator its Initial Value, a state-space block C x₀ + D u₀ — and integrates nothing (ICoreBlockSolverEnvironment::solve, the solvedOnce flag). Only from the second visit on does a continuous block advance its state from the time it was last solved. This is what the exported code does (y = C x + D u from the current state, then the update) and what Simulink reports at t = 0. Until 2026-09-02 the first visit integrated one full period from −h to startTime, so a Constant 1 into an Integrator read 0.1 at t = 0 and startTime + h when the start time was not zero.

The stop time is solved, and nothing beyond it. A finite run ends on the first solve time that is at or within STOP_TIME_SLACK (1e-9 s, ICoreModelSimulator.h) of stopTime; a grid point or sub-step past that is not taken (reachedStopTime, coreLogic, ICoreSimulatorRepeatableWindowSolvingList::fire). So Ts = 0.1, stopTime = 10 solves t = 10 as its last sample in every mode. Before 2026-09-02 the three modes disagreed: the continuous window stopped one grid point short (the collector still wrote a row for the stop time, repeating the previous values — the flat last segment on the older block-page plots), while Discrete and variable step solved one step past it.

Fixed step. Time is a grid, not an accumulator: t_n = startTime + n · h, recomputed from the integer step index every step so that Ts = 0.1 lands the tenth sample on exactly 1.0 and not 0.9999999999999999 (advanceTimeByOneStep, which notes every exported code target computes its clock the same way). Under Discrete, h is Global Sampling Time (sec) (default 0.1). Under Continuous fixed-step the run is multi-rate: every distinct block rate is rounded to the Multi-rate Sampling Tolerance (default 1e-9 s), the greatest common divisor of the rounded rates is the finest sub-step and their least common multiple is one "repeatable window"; each block is solved only on the sub-steps that are multiples of its own rate (initializeSamplingTimes, ICoreSimulatorRepeatableWindow::fire). How a block's rate is resolved is on Sample time and loops — how the simulator paces a diagram.

Implicit fixed step (BE1, TR2; since 2026-09-02). The two methods for stiff models — a fast electrical pole beside a slow mechanical one, a stiff contact or friction law, any block whose fastest mode makes an explicit method's step unusably small. Both are solved per block: each visit finds x_n from its own implicit equation by Newton iteration on G(x) = x − c − γ f(x, u_n, t_n), the Jacobian of f taken by forward differences one state entry at a time, started from an explicit Euler predictor, at most 25 iterations to |δ| ≤ 1e-12 + 1e-10 |x| (solveImplicitStep, ICoreRungeKuttaEstimation.cpp). A linear block converges in two iterations — the first to finite-difference accuracy, the second to rounding — so for a state-space or transfer-function block the answer is the exact implicit step. A step that does not converge keeps its last iterate and is reported once per run ("Implicit step did not converge at block: …"); a smaller Global Sampling Time is the usual fix. Two things to know: Backward Euler with a constant input reproduces the exporter's Backward Euler discretization sample for sample (x_k = (I − hA)⁻¹ (x_{k−1} + hB u)), while for a changing input the simulator uses u_n (the input at the new time) where the exported recursion holds u[k]; and the implicit equation is per block, so stiffness that lives across a feedback loop is still coupled explicitly through the one-sample-stale edge on Sample time and loops — how the simulator paces a diagram. TR2 is the trapezoidal rule (second order, A-stable, and it rings rather than damps on a very stiff mode: x_k = ((1 + hλ/2)/(1 − hλ/2)) x_{k−1} has ratio −49/51 at hλ = −100); BE1 is first order and L-stable (ratio 1/101 there). Prefer BE1 when the fast modes are noise to be killed, TR2 when the slow dynamics must be accurate. Both are fixed-step, so code export accepts them.

Variable step (RK45, RK23). The first step is Initial Time Step (default 0.01). After every step each integrating block proposes the next size from its embedded error estimate (formula below) and the simulator takes the smallest proposal, clamped to [Min Time Step, Max Time Step] (defaults 0.001 and 0.1); the step is then shortened, if need be, to land exactly on the next sample hit of any discrete-only block and exactly on stopTime — calculateNextStepTimeDelta. Time accumulates, t += h. Two things the scheme does not do, on purpose: it never rejects a step — every step is accepted and its error only sizes the next one, because a rejected step would have to roll back every block's state and every port, and blocks with their own memory (scopes, running statistics, random sources) cannot be rolled back generically; and it does not step below Min Time Step for accuracy — when the controller asks for less, the step is clamped, taken, and one warning per run says the tolerance is not met there ("Variable-step solver asked for a step of …"). The two clamps that do go below Min Time Step are the sample hit and the stop time, which are facts about the model rather than preferences about the integrator.

Implicit variable step (TRBDF2; since 2026-09-02) is the stiff counterpart of RK45/RK23 — Simulink's ode23tb. One step is two implicit stages solved per block by the same Newton iteration as BE1/TR2: the trapezoidal rule to t + γh, then second-order BDF to t + h, with γ = 2 − √2, d = γ/2, w = √2/4 (the DIRK tableau c = (0, γ, 1), a = ((0,0,0), (d,d,0), (w,w,d)), b = (w, w, d)). It is L-stable, so a fast mode is damped in one step rather than resolved. The embedded third-order weights b̂ = ((1−w)/3, (3w+1)/3, d/3) give the error e = h Σ (b_i − b̂_i) k_i with k3 = f(x_n, u_n, t_n), and — following Hosea & Shampine — that estimate is filtered through (I − d h J)⁻¹, the matrix the last Newton iteration already formed (filterThroughNewtonMatrix), without which a stiff mode inflates the estimate and keeps the step tiny. The filtered error then goes through the same scaled norm and controller as the explicit pairs, with exponent 1/3; the first step, the sample-hit and stop-time clamps, the Min Time Step warning and the no-rejection rule all apply unchanged. Measured 2026-09-02 on this page's Constant 1 → State Space (A = −1000, B = 1000, C = 1) recipe, initialTimeStep 0.01, minTimeStep 0.001, maxTimeStep 0.1, relativeTolerance 1e-3, absoluteTolerance 1e-6, 0 to 2 s, through docsSample --recipe (build of commit ddc5e94f plus this change):

RK45    610 rows, steps 0.01, 0.0024, then 0.001 (the Min Step Size clamp) for the rest of the run;
        largest step 0.01; y(2) = 0.999561, wandering between 0.997 and 0.9998 after the transient
TRBDF2   28 rows, steps 0.01, 0.002, 0.001, 0.0014, 0.0022, 0.0043, … opening to 0.1;
        y = 1.000000 from t = 0.05 on, y(2) = 1

RK45 is not inaccurate here — it is held at its stability boundary, |hλ| ≈ 1 at the clamp, and pays 600 steps for a solution that is a constant. Code export refuses TRBDF2, as it refuses every variable-step method.

Joint integration — the continuous coupling (since 2026-09-02; the Coupling setting, solverCoupling in the console, default Joint). Under the Continuous solver every block's continuous states are gathered into one state vector and advanced together by the chosen method, with the whole diagram re-evaluated at every intermediate stage (ICoreJointIntegration, Core/ICoreJointIntegration.cpp). One step from t_{n−1} to t_n: the trial state X of a stage is written into every state block, the blocks are swept — first every state block with no direct feedthrough writes the output its trial state alone determines, then the blocks in solve order each write the output computed from their trial state (or start-of-step state when they have none) and their inputs — and then every state block's derivative is read from the ports the sweep left. So the input an Integrator sees at a stage is the output the rest of the diagram produces for that trial state, feedback included, and a loop between continuous blocks is integrated exactly: Integrator → Gain −1 → Integrator under RK4 at h = 0.1 reproduces e^{−t} to 1e-6 where the per-block scheme is off by O(h) (the solver suite). The first phase of the sweep is what makes that true whatever the solve order: the build opens every loop at exactly such a block (Sample time and loops — how the simulator paces a diagram), so its consumers may be ordered before it, and without the phase they read its port one stage stale — the same loop with the Gain created first, or a Scope on the diagram, came out 0.3585 for e⁻¹ = 0.3679 until 2026-09-03. The commit pass does the same: those outputs from the new state first, then everything in solve order. The same holds for the implicit methods, whose Newton iteration then runs over the whole vector — a stiff loop (Integrator → Gain −1000) is implicit as a loop under TRBDF2 and walks its settled solution in a few dozen steps where per-block coupling, explicit through the stale edge, is held at Min Time Step. A source is evaluated at the stage times rather than interpolated between visits (Sine → Integrator gives 1 − cos t to 1e-5 at h = 0.1).

Three things decide who takes part. A block with continuous states is one that says so: hasContinuousStates() is true once initializeStateSpace_Continuous() has run, and a block with a real compute_f but no state space — the parameter-varying filters — calls setHasContinuousStates(true); a block that does not is committed at t_n from its start-of-step state and its derivative is never asked for. A block is swept at the stages only if its output is pure — a function of (x, u, t) with no side effect; one that draws a random number, latches, logs, records or asserts says setOutputHasSideEffects(true) (Backlash, Relay, Playback, Signal Generator, every sink and scope, the model-verification checks) and its ports hold their last committed value through the stages, which keeps such a block on exactly the one-visit-per-step it always had. A discrete-only block never takes part in a sweep; it runs its own recursion at its sample hits in the commit pass. The commit pass, in solve order, is the same single visit as before: outputs at t_n from the new state and the ports as they stand. When it applies: variable step, and single-rate fixed step. The multi-rate window (blocks at different rates) keeps the per-block scheme, because its blocks do not step together. The Discrete solver is untouched. Per-block remains selectable, and it is the scheme every exported code target reproduces — a feedback diagram verified against its export under Joint shows the delay difference the export still has (the ten targets run the single-visit recursion); switch to Per-block for that comparison, or verify under the Discrete solver as the parity suites do.

Measured 2026-09-02 (headless, docsSample --recipe on Integrator (initial value 1) → Gain k → Integrator, i.e. dx/dt = k x; the same binary as the rest of this page):

k = -1,    RK4, Ts = 0.1, 0..1 s        x(1)  Joint: 0.3678798   Per-block: 0.348175   exact e^-1 = 0.3678794
                                          error Joint: 3.3e-7     Per-block: 2.0e-2
k = -1000, TRBDF2, initial 0.01, min 0.001, max 0.1, relTol 1e-3, absTol 1e-6, 0..0.5 s
           Joint:     20 rows, steps 0.01, 0.002, 0.001, ... opening to 0.1; max |x| = 1; x(0.5) = -2e-16
           Per-block:  8 rows, steps 0.01, 0.05, 0.1, 0.1 ...;  max |x| = 6.7e9;  x(0.5) = -6.7e9

The per-block run of the stiff loop is the instructive one: each block's own problem is "integrate a constant", its embedded error is nil, the controller opens the step to the maximum — and the loop, explicit through the stale edge at h · 1000 = 100, diverges without a single warning. No per-block estimate can see a loop.

Known discontinuities under variable step (since 2026-09-02). A source whose output is known in advance to jump or bend announces the next such time through ICoreBlockSolverEnvironment::nextDiscontinuityTimeAfter(t): the Step's step time (one per entry), the Ramp's start time, every breakpoint of a Repeating Sequence (in every period), a Signal Editor or From Workspace table, the delay and every edge of a Waveform Generator's square, pulse or sawtooth, and the edges of a Signal Generator's square or sawtooth. The solver lands the clock on the announced time in two moves — an approach step that ends max(1e-8, 1e-3 · Initial Time Step) s before it (1e-5 s at the default), then a step onto it — and restarts the step size from Initial Time Step afterwards (nextDiscontinuityAfter, eventApproachDelta, calculateNextStepTimeDelta). The approach step is the point: the step that ends on the event is the one whose input interpolation runs from the pre-event to the post-event value, so making it a thousandth of a step makes the smear a thousandth of what it was. Measured in the solver suite: a Step at 1 s into an Integrator under RK45 reads 0 at t = 1 to within 1e-5 and t − 1 after it, where before the landing the integral could carry up to h/2 = 0.05 from the step that happened to contain the jump. What this is not: a zero-crossing locator. A discontinuity that depends on a signal — a Saturation or Relay on a plant output, an Integrator reaching its limit, the Hit Crossing block — is not found inside a step, because locating it means retaking the step, and this simulator does not reject steps (above). Those blocks still switch at whichever step first sees the crossing. Pulse Generator, Repeating Sequence Stair and PWM are discrete-only blocks and land on their sample hits instead.

Linearizing the whole diagram (since 2026-10-01). At the current operating point (x, u), the diagram is linearized to dx/dt = A dx + B du, dy = C dx + D du, as Simulink's linmod does between the root Inports and Outports. Each column of A and C is a central difference over one state, (f(x + h e_j) − f(x − h e_j)) / 2h with h = 1e-5 · max(1, |x_j|), and each column of B and D over one input element, every evaluation being the same sweep and derivative a joint step takes, so a feedback loop is closed exactly as the run closes it. The error is of order h²: on a smooth model about 1e-9, so it agrees with Simulink's analytic, block-by-block Jacobian to a band rather than to the last bit. Only states that move are states: a feed-through block's placeholder is not listed.

Linearizing between points on links (since 2026-10-01). Instead of the root ports, the linearization can be taken between analysis points on links, the points Simulink Control Design places with linio. A point sits on the output port that sources a link and carries three independent flags, exactly as a Simulink output port stores them: input, output and open loop. Wherever the point's block computes its value v during a sweep, three things happen in this order: an output point records v as its element of y; an open-loop point replaces v with the value it had at the operating point, so whatever reads the link no longer closes a loop through it; an input point adds its element of du. Simulink's named types are combinations of the three: inout is input and output, openinput input and open loop, openoutput output and open loop, and loopbreak open loop alone. With input and output on the same point, y is measured before the perturbation. In a recipe a point is set with block.outputPort(i).setLinearization(input, output, openloop), or with one of those names, and none clears it.

Scheduled hits land without the approach step (since 2026-09-30). A block that calls setAnnouncesScheduledHits(true) says the times it announces are hits it schedules — a Hit Scheduler's now + Δt — not jumps in its own output. The solver lands on such a time in one ordinary step and still restarts the step size after it. That is Simulink's behaviour, measured: under ode45 a hit requested at 1.2345 s is reached straight from 1.0. A time that any unflagged block also announces keeps its approach step. The question is asked after every step is solved, so a block may compute its next hit from the input it read in compute_h and keep it in its own state. Pinned by the solver case variable_step_lands_on_a_scheduled_hit_without_an_approach. Code export refuses such a block by name, because an exported core runs a fixed step and cannot land on a time computed during the run. The library's block of this kind is Control_Systems/Messages_And_Events/Hit_Scheduler. Simulink's block inherits its enable's rate, so it reads the enable and schedules a hit only on that rate's sample hits. ICore applies the same rule when the enable comes straight from a discrete-only block: a Repeating Sequence Stair at 0.1 s, true at 0, 0.1 and 0.2, with Δt = 0.35, gives hits at 0, 0.35, 0.45 and 0.55 under ode45 in R2026a and under RK45 in ICore. An enable from any other block is read on every solved step. Pinned by run_analysis_blocks.

Discrete-only blocks under variable step (unit delays, discrete filters, zero-order holds — every block that declares setDiscreteOnlyBlock(true)) keep their own period: they are solved only when the clock is on one of their sample hits, startTime + k · Ts within STOP_TIME_SLACK, and hold their ports in between (fireGlobalSolve, isOnSampleHit, nextSampleHitAfter). Before 2026-09-02 they were stepped at whatever size the integrator picked. Continuous blocks with a positive Sampling Time (s) are still stepped every variable step, with the warning on Sample time and loops — how the simulator paces a diagram.

Infinite runs. With infiniteSimulation on, the clock is rewound to startTime once it has run more than maximumConsecutiveTimeBuffer (default 10000 s) past it. The rewind moves every block's own clock by the same amount (shiftLastSolvedTime), so the next step is one ordinary period; states are kept — a rewind is a clock event, not a reset, and time-based sources (Step, Sine, Clock) start their timeline again. Until 2026-09-02 only the global clock was rewound, and every continuous block was handed a step of minus the whole span.

The integrators (Core/ICoreRungeKuttaEstimation.cpp). For a continuous block, f(x, u, t) is the block's state derivative — for a state-space block literally A x + Bu u (+ Bf f) — and one visit advances x from t_{n-1} to t_n = t_{n-1} + h. The block's inputs at the two ends, u_{n-1} and u_n, are known; inputs at intermediate stages are linearly interpolated between them, u(c) = (1 − c) u_{n-1} + c u_n, and the output is then y_n = h(x_n, u_n, t_n) from the updated state and the current input. With k1 = f(x_{n-1}, u_{n-1}, t_{n-1}):

RK1  (Euler)          x_n = x_{n-1} + h k1
RK2  (Heun)           k2 = f(x_{n-1} + h k1, u_n, t_n)
                      x_n = x_{n-1} + (h/2)(k1 + k2)
RK3                   k2 = f(x_{n-1} + (h/2) k1,          u(1/2), t + h/2)
                      k3 = f(x_{n-1} − h k1 + 2h k2,     u_n,    t_n)
                      x_n = x_{n-1} + (h/6)(k1 + 4 k2 + k3)        (Kutta's third-order rule)
RK4  (classical)      k2 = f(x_{n-1} + (h/2) k1, u(1/2), t + h/2)
                      k3 = f(x_{n-1} + (h/2) k2, u(1/2), t + h/2)
                      k4 = f(x_{n-1} + h k3,     u_n,    t_n)
                      x_n = x_{n-1} + (h/6)(k1 + 2 k2 + 2 k3 + k4)
BE1  (Backward Euler)  x_n = x_{n-1} + h f(x_n, u_n, t_n)                      (implicit; Newton per block)
TR2  (Trapezoidal)     x_n = x_{n-1} + (h/2) [ k1 + f(x_n, u_n, t_n) ]         (implicit; Newton per block)
RK23 (Bogacki–Shampine)  stages at c = 1/2, 3/4, 1 with the standard tableau;
                      x_n = third-order solution (2/9, 1/3, 4/9);
                      e   = x_n − second-order solution (7/24, 1/4, 1/3, 1/8)
                      err = RMS_i ( e_i / (absTol + relTol · max(|x_{n-1,i}|, |x_{n,i}|)) )
                      h_next = h · clamp( 0.9 · err^(−1/3), 0.2, 5 )
RK45 (Dormand–Prince) seven stages at c = 1/5, 3/10, 4/5, 8/9, 1, 1 with the standard tableau;
                      x_n = fifth-order solution;  e = x_5th − x_4th
                      err = RMS_i ( e_i / (absTol + relTol · max(|x_{n-1,i}|, |x_{n,i}|)) )
                      h_next = h · clamp( 0.9 · err^(−1/5), 0.2, 5 )
TRBDF2 (implicit)     γ = 2 − √2, d = γ/2, w = √2/4
                      x_γ = x_{n-1} + d h [ k1 + f(x_γ, u(γ), t + γh) ]              (trapezoidal, Newton)
                      x_n = x_{n-1} + w h (k1 + k2) + d h f(x_n, u_n, t_n)          (BDF2, Newton)
                      e   = (I − d h J)⁻¹ · h [ ((4w−1)/3) k1 − (1/3) k2 + (2d/3) k3 ]
                      err as above;  h_next = h · clamp( 0.9 · err^(−1/3), 0.2, 5 )

relTol and absTol are Relative Tolerance (default 1e-3) and Absolute Tolerance (default 1e-12), with their standard meanings: the error of each state entry is measured against absTol + relTol · |x|, so err ≤ 1 means every entry is inside its band, and the next step is the current one scaled by 0.9 · err^(−1/(q+1)) — q the order of the lower solution — bounded to between one fifth and five times the current step (scaledErrorNorm, nextStepFactor). A zero error grows the step by the full factor of five; a NaN shrinks it by five. (Until 2026-09-02 the norm was unscaled — relTol acted as an absolute error target and absTol as a floor under it — and the growth per step was unbounded below Max Time Step.) RK3 computes Kutta's classical third-order rule; the genuine Bogacki–Shampine pair is RK23. An unrecognised stepping string stops the run with "Invalid solver type".

Where the one-sample delay comes from. A block's inputs are read from its source ports as they stand when the block is visited, in solve order, so a signal that comes from later in the order (a feedback edge) is the previous step's value. That is a property of the commit pass, not of any integrator, and it is described once, with a measured Simulink comparison, on Sample time and loops — how the simulator paces a diagram. Under Joint coupling it no longer affects the continuous states — their derivatives are read after a full sweep at every stage — but it still governs a loop closed through a discrete or side-effecting block, and the Discrete solver.

What step, impulse and ramp compute#

ICoreTimeResponse::compute() (src/ICoreBlocks/ICoreMath/ControlSystems/StaticMethods/ICoreTimeResponse.cpp) backs the console's own three-argument forms step(G, duration, points) and impulse(G, duration, points), ramp(G[,duration[,points]]), and the member forms G.step(...), G.impulse(...), G.ramp(...). (step(sys), step(sys, Tfinal), step(sys, t) and the same forms of impulse are MATLAB's: they answer the output y alone, as a column, on MATLAB's own time grid, and [y, t] = step(...) returns the grid too — see Command glossary — console commands, verbs, functions.) It is always a discrete simulation from zero initial state:

1. G is realized as a state space by fromTransferFunction() (the canonical form above).
2. If G is continuous, the grid is Ts = duration / (points − 1) and the realization is ZOH-discretized at that Ts — always ZOH, whatever the solver setting says. If G is discrete, the grid is its own Ts, points is ignored and the count is floor(duration / Ts) + 1 (at least 2, at most 2 000 000).
3. Then, for k = 0 … N−1, t_k = k Ts; y_k = Cd x_k + Dd u_k; x_{k+1} = Ad x_k + Bd u_k with u_k = 1 (step), u_k = t_k (ramp), or the impulse u_0 = 1/Ts for a continuous G — a unit sample of width Ts has area Ts, so it is scaled to unit area — and u_0 = 1 for a discrete one, u_k = 0 afterwards.

The result is an N×2 [time, output] time series. Because the output uses x_k before the update, y_0 = Dd u_0 — zero for a strictly proper G — and the first non-zero sample is at t = Ts. When no duration is given, suggestedDuration() picks one from the poles: five time constants of the slowest decaying mode (5 / min|Re s|, with a discrete pole mapped through s = ln(z)/Ts), or three periods of the slowest oscillation when nothing decays, or 10 s when the poles say nothing; the value is rounded up to the next 1/2/5 × 10^k, and for a discrete model clamped to between 20 and 20 000 samples. The default point count is 500 (kDefaultResponsePoints, ICoreExpressionEvaluator.cpp), which is why G.step() on 1/(s+2) below answers with 500 samples.

Real run — a first-order system, by hand and by the console#

Commands run 2026-09-26 against the built app (binary of 2026-09-25 23:58, approximately commit 1dac24c57), one process per line as HOME=<scratch> ICoreBlocks.app/Contents/MacOS/ICoreBlocks --console "<line>"; startup lines stripped, everything else verbatim. The system is G = 1/(s+2), so a = −2, Ts = 0.1.

>> G = tf([1],[1 2]); G
G = 
    1
  -----
  s + 2
Continuous-time transfer function.  # Transfer Function
>> G = tf([1],[1 2]); tf2ss(G)
  A =
       x1
   x1  -2
  B =
       u1
   x1   1
  C =
       x1
   y1   1
  D =
       u1
   y1   0
Continuous-time state-space model.  # State Space
>> H = tf([1 3],[1 2 5]); tf2ss(H)
  A =
       x1  x2
   x1  -2  -5
   x2   1   0
  B =
       u1
   x1   1
   x2   0
  C =
       x1  x2
   y1   1   3
  D =
       u1
   y1   0
Continuous-time state-space model.  # State Space

The ZOH matrices, computed the way the discretizer computes them — the augmented exponential — and the scalar closed form beside it:

>> expm([-2 1; 0 0]*0.1)
[[0.818731, 0.0906346], [0, 1]]  # Matrix of Double
>> expm(-2*0.1)
0.818731  # Double

So Ad = e^{−0.2} = 0.818731 and Bd = (Ad − 1)/a = 0.0906346. The step and impulse responses over 1 s in 11 points (Ts = 0.1):

>> G = tf([1],[1 2]); Y = step(G, 1, 11); Y.time()
[[0], [0.1], [0.2], [0.3], [0.4], [0.5], [0.6], [0.7], [0.8], [0.9], [1]]  # Matrix of Double
>> G = tf([1],[1 2]); Y = step(G, 1, 11); Y.values()
[[0], [0.0906346], [0.16484], [0.225594], [0.275336], [0.31606], [0.349403], [0.376702], [0.399052], [0.417351], [0.432332]]  # Matrix of Double
>> G = tf([1],[1 2]); Y = impulse(G, 1, 11); Y.values()
[[0], [0.906346], [0.742054], [0.607542], [0.497413], [0.407248], [0.333426], [0.272986], [0.223502], [0.182988], [0.149818]]  # Matrix of Double
>> G = tf([1],[1 2]); G.step()
Step response plotted (500 samples)

y_1 = Bd = 0.0906346 and y_2 = Ad·y_1 + Bd = 0.16484 are the recursion of step 3 above; the impulse's first sample is Bd · (1/Ts) = 0.906346, then it decays by Ad each sample (0.742054 = 0.818731 × 0.906346). Cross-checked in Python (stdlib only, same date):

$ python3 -c '
import math
a,b,Ts=-2.0,1.0,0.1
Ad=math.exp(a*Ts); Bd=(Ad-1)/a*b
print("Ad =",round(Ad,6),"Bd =",round(Bd,7))
x=0.0; ys=[]
for k in range(11): ys.append(x); x=Ad*x+Bd
print("step:",[round(v,6) for v in ys])
x=0.0; yi=[]
for k in range(11): yi.append(x); x=Ad*x+Bd*(1/Ts if k==0 else 0.0)
print("impulse:",[round(v,6) for v in yi])
print("tustin Ad =",round((1+a*Ts/2)/(1-a*Ts/2),6),"fwd =",round(1+a*Ts,6),"bwd =",round(1/(1-a*Ts),6))'
Ad = 0.818731 Bd = 0.0906346
step: [0.0, 0.090635, 0.16484, 0.225594, 0.275336, 0.31606, 0.349403, 0.376702, 0.399052, 0.417351, 0.432332]
impulse: [0.0, 0.906346, 0.742054, 0.607542, 0.497413, 0.407248, 0.333426, 0.272986, 0.223502, 0.182988, 0.149818]
tustin Ad = 0.818182 fwd = 0.8 bwd = 0.833333

Every printed sample agrees with the console to the six digits it prints. The last line is the same pole under three other methods from the table — Tustin 0.818182, Forward Euler 0.8, Backward Euler 0.833333 against the exact e^{−0.2} = 0.818731 — which is the size of the choice at Ts = one fifth of the time constant. The solver settings the run used, read back from the same binary:

>> getModelConfig solverType
Continuous
>> getModelConfig steppingType
Fixed-step | Runge-Kutta | RK4
>> getModelConfig discretizationMethod
Zero-order Hold
>> getModelConfig globalSamplingTime
0.1

For contributors#

The contributor pages on the numerics module, the simulator and the block model (icoremath, icoresimulation, icoremodel in the front-matter related: list) hold the invariants this mathematics lives under; the console regression corpus that pins these numbers is testingLabs/tests/regression/controlsystems.itest, and the stepping scheme itself — first visit, stop time, variable-step first step, growth bound and sample hits — is pinned by the solver suite (testingLabs/RegressionLab/ICoreRegressionSuite_Solver.cpp, regressionTest --suite solver).