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/ICoreSDK/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/Continues, Discrete State Space — Control Systems/Discrete — are the special case f = 0. The Nonlinear State Space — Control Systems/Continues 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 (isContinues()); 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 poly([..]) 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()), and the numerator's order must not exceed it — a non-causal quotient is refused with "Invalid numerator and/or denominator orders" (tryUpdatingCoefficients). 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. A zero (or < 1e-15) leading denominator coefficient, or an order of 0, is refused ("Invalid transfer function") 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/ICoreSDK/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/ICoreSDK/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/ICoreSDK/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 | Bogacki–Shampine | RK3, Fixed-step | Runge-Kutta | RK4 (the default), Variable-step | Dormand-Prince | RK45, Variable-step | Bogacki–Shampine | RK23. A Discrete solver is always fixed-step.

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.

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 and the simulator takes the smallest proposal, clamped to [Min Time Step, Max Time Step] (defaults 0.001 and 0.1) — calculateNextStepTimeDelta. Time then accumulates, t += h.

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)
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);
                      err = ‖x_n − second-order solution (7/24, 1/4, 1/3, 1/8)‖_RMS
                      h_next = h · (relTol / max(err, absTol))^(1/3)
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;  err = ‖x_5th − x_4th‖_RMS
                      h_next = h · (relTol / max(err, absTol))^(1/5)

relTol and absTol are Relative Tolerance (default 1e-3) and Absolute Tolerance (default 1e-12); the error norm is the RMS over the state vector's entries. Note that the RK3 entry is labelled "Bogacki–Shampine" in the setting's name but 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 scheme, 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.

What step, impulse and ramp compute#

ICoreTimeResponse::compute() (src/ICoreSDK/ICoreMath/ControlSystems/StaticMethods/ICoreTimeResponse.cpp) backs the console's step(G[,duration[,points]]), impulse(...), ramp(...) and their G.step(...) forms. 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 step(G) on 1/(s+2) below answers with 500 samples over 0..5 s.

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

Commands run 2026-08-17 against the built app (binary of 2026-08-16 23:39, approximately commit 2e126fbf), one process per line as HOME=<scratch> QT_QPA_PLATFORM=offscreen build-mac/ICoreBlocks.app/Contents/MacOS/ICoreBlocks --console "<line>"; startup noise 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)  # Transfer Function
>> G = tf([1],[1 2]); tf2ss(G)
A = [-2]
B = [1]
C = [1]
D = [0]
  # State Space
>> H = tf([1 3],[1 2 5]); tf2ss(H)
A = [-2 -5; 1 0]
B = [1; 0]
C = [1 3]
D = [0]
  # 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]); step(G)
Time Series (500 samples, t = 0 .. 5)  # Time Series

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.