function [y, t_ret, x_arr] = __linear_simulation__ (sys, u, t, x0) method = "zoh"; [urows, ucols] = size (u); if (isct (sys)) # continuous-time system if (isempty (t)) # lsim (sys, u, [], ...) error ("lsim: invalid time vector"); else len_t = length (t); if (len_t == 1) # lsim (sys, u, tfinal, ...) len_t = urows; tsam = t / (len_t - 1); tinitial = 0; tfinal = t; else # lsim (sys, u, t, ...) if (urows != len_t) error ("lsim: input vector u must have %d rows", len_t); endif tinitial = t(1); tfinal = t(end); t_ret = linspace (tinitial,tfinal,len_t)' ; dt_diff = t_ret - reshape (t,[],1); # check if t is regularly spaced if ( max(abs(dt_diff)) > 1e3*eps ) # tolerance strategy may be disputed ! error ("lsim: time vector t must be regularly spaced"); endif tsam = (t(end) - t(1)) / (len_t-1) ; # safe : t has been checked to be regularly spaced endif endif tsam_NS = .5 * pi / max (abs (eig (sys))); # generalized Nyquist-Shannon multi_NS = .5; # may be disputed ! if (tsam >= multi_NS * tsam_NS) warning ("lsim: sampling of continuous-time system is too coarse"); endif sys = c2d (sys, tsam, method); # convert to discrete-time model else # discrete-time system len_t = urows; tsam = abs (get (sys, "tsam")); # use 1 time unit as default if tsam is unspecified (-1) if (isempty (t)) # lsim (sys, u) toff = 0; else # lsim (sys, u, toffset) if (1 != length (t)) error ("lsim: discrete-time systems only allow for a scalar time to be supplied, namely the offset"); endif toff = t; endif tinitial = toff; tfinal = toff + tsam * (len_t-1) ; endif [A, B, C, D] = ssdata (sys); [p, m] = size (D); # number of outputs and inputs n = rows (A); # number of states ## (reconstructed) time vector to return if (1 != exist ("t_ret", "var") ) # the work may have been done already t_ret = linspace (tinitial,tfinal,len_t)' ; endif if (ucols != m) error ("lsim: input vector u must have %d columns", m); endif ## preallocate memory y = zeros (len_t, p); x_arr = zeros (len_t, n); ## initial conditions if (isempty (x0)) x0 = zeros (n, 1); elseif (n != length (x0) || ! is_real_vector (x0)) error ("lsim: x0 must be a vector with %d elements", n); endif x = reshape (x0, [], 1); # make sure that x is a column vector ## simulation for k = 1 : len_t y(k, :) = C * x + D * u(k, :).'; x_arr(k, :) = x; x = A * x + B * u(k, :).'; endfor endfunction