Coverage for src / sdynpy / core / sdynpy_nonlinear_system.py: 9%
332 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-11 16:22 +0000
1# -*- coding: utf-8 -*-
2"""
3Defines a nonlinear system using a generic equation of motion.
5The SystemNL consists of a generic equation of motion that defines its behavior.
6The SystemNL also contains a CoordinateArray defining the degrees of freedom of
7the System.
8"""
9"""
10Copyright 2022 National Technology & Engineering Solutions of Sandia,
11LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S.
12Government retains certain rights in this software.
14This program is free software: you can redistribute it and/or modify
15it under the terms of the GNU General Public License as published by
16the Free Software Foundation, either version 3 of the License, or
17(at your option) any later version.
19This program is distributed in the hope that it will be useful,
20but WITHOUT ANY WARRANTY; without even the implied warranty of
21MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22GNU General Public License for more details.
24You should have received a copy of the GNU General Public License
25along with this program. If not, see <https://www.gnu.org/licenses/>.
26"""
28import numpy as np
29from scipy.interpolate import interp1d
30from scipy.integrate import solve_ivp
31from scipy.optimize import approx_fprime
32from scipy.linalg import eig
33from scipy.signal import butter, filtfilt
34from .sdynpy_coordinate import CoordinateArray,coordinate_array
35from .sdynpy_system import System
36from .sdynpy_data import TimeHistoryArray,time_history_array
37from .sdynpy_shape import shape_array
38from ..signal_processing import generator
40class SystemNL:
41 """Function representing a Structural Dynamics System"""
43 def __init__(self, coordinate : CoordinateArray,
44 state_derivative_function,
45 state_size = None,
46 displacement_output_function = None,
47 velocity_output_function = None,
48 acceleration_output_function = None,
49 excitation_output_function = None):
50 """Creates a nonlinear system object
52 Parameters
53 ----------
54 coordinate : CoordinateArray
55 The physical degrees of freedom associated with the system
56 state_derivative_function : function
57 A function with calling signature (time, state, input) that returns
58 the state derivative.
59 state_size : int, optional
60 The number of state degrees of freedom in the system. If not provided,
61 it will be assumed to be the same size as coordinate.
62 displacement_output_function : _type_, optional
63 _description_, by default None
64 velocity_output_function : _type_, optional
65 _description_, by default None
66 acceleration_output_function : _type_, optional
67 _description_, by default None
68 excitation_output_function : _type_, optional
69 _description_, by default None
70 """
71 self.coordinate = coordinate.flatten()
72 if state_size is None:
73 self.state_size = self.coordinate.size*2
74 else:
75 self.state_size = state_size
76 self.state_derivative_function = state_derivative_function
77 if displacement_output_function is None:
78 def default_disp_fn(sys_state, sys_input, response_coordinate = None):
79 if response_coordinate is not None:
80 indices,signs = self.coordinate.find_indices(response_coordinate)
81 else:
82 indices = slice(None)
83 signs = np.array([1])
84 return signs[:,np.newaxis]*sys_state[:self.state_size//2][indices]
85 self.displacement_output_function = default_disp_fn
86 else:
87 self.displacement_output_function = displacement_output_function
88 if velocity_output_function is None:
89 def default_vel_fn(sys_state, sys_input, response_coordinate = None):
90 if response_coordinate is not None:
91 indices,signs = self.coordinate.find_indices(response_coordinate)
92 else:
93 indices = slice(None)
94 signs = np.array([1])
95 return signs[:,np.newaxis]*sys_state[self.state_size//2:][indices]
96 self.velocity_output_function = default_vel_fn
97 else:
98 self.velocity_output_function = velocity_output_function
99 if acceleration_output_function is None:
100 def default_accel_fn(time, sys_state, sys_input, response_coordinate = None):
101 if response_coordinate is not None:
102 indices,signs = self.coordinate.find_indices(response_coordinate)
103 else:
104 indices = slice(None)
105 signs = np.array([1])
106 return signs[:,np.newaxis]*self.state_derivative_function(time, sys_state,sys_input)[self.state_size//2:][indices]
107 self.acceleration_output_function = default_accel_fn
108 else:
109 self.acceleration_output_function = acceleration_output_function
110 if excitation_output_function is None:
111 def default_ex_fn(sys_state, sys_input, excitation_coordinate = None):
112 if excitation_coordinate is not None:
113 indices,signs = self.coordinate.find_indices(excitation_coordinate)
114 else:
115 indices = slice(None)
116 signs = np.array([1])
117 return signs[:,np.newaxis]*sys_input[indices]
118 self.excitation_output_function = default_ex_fn
119 else:
120 self.excitation_output_function = excitation_output_function
122 @property
123 def num_dof(self):
124 return self.coordinate.size
126 @classmethod
127 def from_linear_system(cls,system : System):
128 if system.massless_dofs.sum() > 0:
129 raise NotImplementedError('Linear systems with massless degrees of freedom are not yet implemented!')
130 A,B,C,D = system.to_state_space()
131 coordinates = system.coordinate
132 nc = coordinates.shape[0]
133 def state_derivative_function(time,sys_state,sys_input):
134 return A@sys_state + B@sys_input
135 def displacement_output_function(sys_state,sys_input, response_coordinate = None):
136 if response_coordinate is not None:
137 indices,signs = coordinates.find_indices(response_coordinate)
138 else:
139 indices = slice(None)
140 signs = np.array([1])
141 return (signs[:,np.newaxis]*C[0*nc:1*nc,:][indices])@sys_state + (signs[:,np.newaxis]*D[0*nc:1*nc,:][indices])@sys_input
142 def velocity_output_function(sys_state,sys_input, response_coordinate = None):
143 if response_coordinate is not None:
144 indices,signs = coordinates.find_indices(response_coordinate)
145 else:
146 indices = slice(None)
147 signs = np.array([1])
148 return (signs[:,np.newaxis]*C[1*nc:2*nc,:][indices])@sys_state + (signs[:,np.newaxis]*D[1*nc:2*nc,:][indices])@sys_input
149 def acceleration_output_function(time,sys_state,sys_input, response_coordinate = None):
150 if response_coordinate is not None:
151 indices,signs = coordinates.find_indices(response_coordinate)
152 else:
153 indices = slice(None)
154 signs = np.array([1])
155 return (signs[:,np.newaxis]*C[2*nc:3*nc,:][indices])@sys_state + (signs[:,np.newaxis]*D[2*nc:3*nc,:][indices])@sys_input
156 def excitation_output_function(sys_state,sys_input, excitation_coordinate = None):
157 if excitation_coordinate is not None:
158 indices,signs = coordinates.find_indices(excitation_coordinate)
159 else:
160 indices = slice(None)
161 signs = np.array([1])
162 return (signs[:,np.newaxis]*C[3*nc:4*nc,:][indices])@sys_state + (signs[:,np.newaxis]*D[3*nc:4*nc,:][indices])@sys_input
163 return cls(coordinates, state_derivative_function,
164 A.shape[0], displacement_output_function,
165 velocity_output_function, acceleration_output_function,
166 excitation_output_function)
168 @classmethod
169 def duffing_oscillator(cls, delta, alpha, beta, coordinate = None):
170 if coordinate is None:
171 coordinate = coordinate_array([1],'X+')
172 def state_derivative_function(time, sys_state, sys_input):
173 return np.array([
174 sys_state[1],
175 - delta*sys_state[1] - alpha*sys_state[0] - beta*sys_state[0]**3 + sys_input[0]
176 ])
177 return cls(coordinate, state_derivative_function)
179 @classmethod
180 def elastic_pendulum(cls, mass, stiffness, unstretched_length, gravity=9.81, node = None):
181 if node is None:
182 node = 1
183 coordinate = coordinate_array(node,['X+','Y+'])
184 m = mass
185 k = stiffness
186 l_0 = unstretched_length
187 g = gravity
188 def state_derivative_function(time, sys_state, sys_input):
189 # State vector x,y,xd,yd
190 x = sys_state[0]
191 y = sys_state[1]
192 xd = sys_state[2]
193 yd = sys_state[3]
194 xdd = k*l_0*x/(m*np.sqrt((l_0 - y)**2 + x**2)) - k*x/m + sys_input[0]
195 ydd = -g - k*l_0**2/(m*np.sqrt((l_0 - y)**2 + x**2)) + k*l_0/m + k*l_0*y/(m*np.sqrt((l_0 - y)**2 + x**2)) - k*y/m + sys_input[1]
196 # Derivative vector xd, yd, xdd, ydd
197 return np.array([xd,yd,xdd,ydd])
198 return cls(coordinate, state_derivative_function)
200 @classmethod
201 def polynomial_stiffness_damping(cls, coordinate, mass, stiffnesses, dampings=None, transformation=None):
202 if dampings is None:
203 dampings = ()
204 num_states = mass.shape[0]*2
205 A_dd = np.concatenate(
206 [-np.linalg.solve(mass,stiffness) for exponent,stiffness in stiffnesses.items()]
207 + ([-np.linalg.solve(mass,damping) for exponent,damping in dampings.items()] if dampings is not None else []),axis=-1)
208 B_dd = np.linalg.solve(mass,transformation.T)
209 def state_derivative(t,sys_state,sys_input):
210 # Assemble the state derivative from the different pieces
211 x = sys_state[:sys_state.shape[0]//2]
212 xd = sys_state[sys_state.shape[0]//2:]
213 z = np.concatenate(
214 [x**exponent for exponent in stiffnesses]
215 + ([xd**exponent for exponent in dampings] if dampings is not None else []))
216 xdd = A_dd@z+B_dd@sys_input
217 return np.concatenate((xd,xdd))
218 def displacement_output_function(sys_state,sys_input, response_coordinate = None):
219 if response_coordinate is not None:
220 indices,signs = coordinate.find_indices(response_coordinate)
221 else:
222 indices = slice(None)
223 signs = np.array([1])
224 return (signs[:,np.newaxis]*transformation[indices])@sys_state[:sys_state.shape[0]//2]
225 def velocity_output_function(sys_state,sys_input, response_coordinate = None):
226 if response_coordinate is not None:
227 indices,signs = coordinate.find_indices(response_coordinate)
228 else:
229 indices = slice(None)
230 signs = np.array([1])
231 return (signs[:,np.newaxis]*transformation[indices])@sys_state[sys_state.shape[0]//2:]
232 def acceleration_output_function(time,sys_state,sys_input, response_coordinate = None):
233 if response_coordinate is not None:
234 indices,signs = coordinate.find_indices(response_coordinate)
235 else:
236 indices = slice(None)
237 signs = np.array([1])
238 dstate = state_derivative(time,sys_state,sys_input)
239 return (signs[:,np.newaxis]*transformation[indices])@dstate[dstate.shape[0]//2:]
240 def excitation_output_function(sys_state,sys_input, excitation_coordinate = None):
241 if excitation_coordinate is not None:
242 indices,signs = coordinate.find_indices(excitation_coordinate)
243 else:
244 indices = slice(None)
245 signs = np.array([1])
246 return (signs[:,np.newaxis]*np.eye(sys_input.shape[0])[indices])@sys_input
247 return cls(coordinate,state_derivative,num_states,displacement_output_function,
248 velocity_output_function,acceleration_output_function,
249 excitation_output_function)
252 def time_integrate(self, forces, dt = None, responses = None, references = None,
253 displacement_derivative = 2, initial_state = None, method = 'RK45',
254 force_interpolation='linear', output_times = None, return_ivp_result = False, **solve_ivp_kwargs):
255 # Set up the inputs and the functions
256 if isinstance(forces, TimeHistoryArray):
257 forces = forces.reshape(-1)
258 dt = forces.abscissa_spacing
259 references = forces.coordinate[...,0]
260 forces = forces.ordinate
261 else:
262 forces = np.atleast_2d(forces)
263 if dt is None:
264 raise ValueError('`dt` must be specified if `forces` is not a `TimeHistoryArray`')
265 # Construct the entire force vector over time
266 if references is None:
267 force_array = forces
268 else:
269 force_array = np.zeros((self.coordinate.shape[0],forces.shape[-1]))
270 reference_indices, reference_signs = self.coordinate.find_indices(references)
271 force_array[reference_indices] = reference_signs[:,np.newaxis]*forces
272 # Create interpolator
273 force_abscissa = np.arange(forces.shape[-1])*dt
274 if output_times is None:
275 output_times = force_abscissa
276 force_interpolator = interp1d(force_abscissa,
277 force_array,kind=force_interpolation,
278 assume_sorted=True, fill_value=(force_array[:,0],force_array[:,-1]),
279 bounds_error=False)
280 def integration_function(t,y):
281 return self.state_derivative_function(t, y, force_interpolator(t))
282 if initial_state is None:
283 initial_state = np.zeros(self.state_size)
284 result = solve_ivp(integration_function, (np.min(output_times),np.max(output_times)),
285 initial_state, method, output_times, **solve_ivp_kwargs)
286 # Package the forces into a time history array
287 forces_ordinate = force_interpolator(output_times)
288 if not result.success:
289 print('Integration was not successful: {:}'.format(result.message))
290 output_times = output_times[...,:result.y.shape[-1]]
291 forces_ordinate = forces_ordinate[...,:result.y.shape[-1]]
292 excitation_ordinate = self.excitation_output_function(result.y,forces_ordinate,references)
293 excitation_signals = time_history_array(
294 output_times,
295 excitation_ordinate,
296 self.coordinate[:,np.newaxis] if references is None else references[:,np.newaxis])
297 if displacement_derivative == 0:
298 fn = self.displacement_output_function
299 elif displacement_derivative == 1:
300 fn = self.velocity_output_function
301 elif displacement_derivative == 2:
302 fn = lambda state,force,response : self.acceleration_output_function(output_times,state,force,response)
303 else:
304 raise ValueError('Displacement Derivative must be one of 0, 1, or 2.')
305 response_ordinate = fn(result.y,forces_ordinate,responses)
306 response_signals = time_history_array(
307 output_times,
308 response_ordinate,
309 self.coordinate[:,np.newaxis] if responses is None else responses[:,np.newaxis])
310 return_vals = [response_signals,excitation_signals]
311 if return_ivp_result:
312 return_vals.append(result)
313 return tuple(return_vals)
315 def __repr__(self):
316 return 'SystemNL with {:} Coordinates ({:} internal states)'.format(self.num_dof, self.state_size)
318 def local_stiffness(self, state_point = None, delta_state = np.float64(1.4901161193847656e-08), time=0):
319 if state_point is None:
320 state_point = np.zeros(self.state_size)
321 forces = np.zeros(self.coordinate.size)
322 fn = lambda state: self.state_derivative_function(time,state,forces)
323 return -approx_fprime(state_point,fn,delta_state)
325 def eigensolution(self, state_point = None, delta_state = np.float64(1.4901161193847656e-08), time=0, group_conjugate_pairs = True):
326 stiffness = self.local_stiffness(state_point, delta_state, time)
327 evals, evects = eig(-stiffness)
328 try:
329 if group_conjugate_pairs:
330 # print('Pairing Complex Conjugate Modes...')
331 evals_reduced = []
332 evects_reduced = []
333 evals_handled = []
334 for index,(evl,evc) in enumerate(zip(evals,evects.T)):
335 if index in evals_handled:
336 continue
337 matches = np.where(np.isclose(evals,evl.conj()) & (np.all(np.isclose(evects.T,evc.conj()),axis=-1) | np.all(np.isclose(-evects.T,evc.conj()),axis=-1)))[0]
338 matches = [match for match in matches if match!=index]
339 if len(matches) > 1:
340 raise ValueError(f'More than one conjugate eigenvalue was found for index {index}: {matches}')
341 if len(matches) == 0:
342 raise ValueError(f'No conjugate eigenvalue was found for index {index}')
343 # else:
344 # print(f' Eigenvalue {index} {evl} matches {matches[0]}')
345 # Choose which to keep
346 if evl.imag > evals[matches[0]].imag:
347 evals_reduced.append(evl)
348 evects_reduced.append(evc)
349 else:
350 evals_reduced.append(evals[matches[0]])
351 evects_reduced.append(evects[:,matches[0]])
352 evals_handled.append(index)
353 evals_handled.append(matches[0])
354 evals = np.array(evals_reduced)
355 evects = np.array(evects_reduced).T
356 except ValueError:
357 print('Could not group conjugate pairs of modes... continuing.')
358 # Compute natural frequency, damping, and mode shape
359 angular_frequencies = np.abs(evals)
360 frequencies = angular_frequencies/(2*np.pi)
361 angular_frequencies[angular_frequencies==0] = 1
362 dampings = -evals.real/angular_frequencies
363 shape_matrix = self.displacement_output_function(evects,np.zeros((self.coordinate.size,evects.shape[-1]))).T
364 shapes = shape_array(self.coordinate,shape_matrix,frequencies,dampings)
365 shapes.sort()
366 return shapes
368 def simulate_test(
369 self, # The system itself
370 bandwidth,
371 frame_length,
372 num_averages,
373 excitation,
374 references,
375 responses=None, # All Responses
376 excitation_level=1.0,
377 excitation_noise_level=0.0,
378 response_noise_level=0.0,
379 steady_state_time=0.0,
380 excitation_min_frequency=None,
381 excitation_max_frequency=None,
382 signal_fraction=0.5,
383 extra_time_between_frames=0.0,
384 integration_oversample=10,
385 displacement_derivative=2,
386 antialias_filter_cutoff_factor=3,
387 antialias_filter_order=4,
388 multihammer_impact_spacing_factor = 4,
389 time_integrate_kwargs = {},
390 generator_kwargs = {}
391 ):
392 available_excitations = ['pseudorandom', 'random',
393 'burst random', 'chirp', 'hammer', 'multi-hammer', 'sine']
394 if not excitation.lower() in available_excitations:
395 raise ValueError('Excitation must be one of {:}'.format(available_excitations))
396 # Create the input signal
397 num_signals = references.size
398 sample_rate = bandwidth * 2 * integration_oversample
399 dt = 1 / sample_rate
400 frame_time = dt * frame_length * integration_oversample
401 df = 1 / frame_time
402 # Create the signals
403 if excitation.lower() == 'pseudorandom':
404 if num_signals > 1:
405 print('Warning: Pseudorandom generally not recommended for multi-reference excitation.')
406 kwargs = {'fft_lines': frame_length // 2,
407 'f_nyq': bandwidth,
408 'signal_rms': excitation_level,
409 'min_freq': excitation_min_frequency,
410 'max_freq': excitation_max_frequency,
411 'integration_oversample': integration_oversample,
412 'averages': num_averages + int(np.ceil(steady_state_time / frame_time))}
413 kwargs.update(generator_kwargs)
414 signals = np.array([
415 generator.pseudorandom(**kwargs)[1]
416 for i in range(num_signals)
417 ])
418 elif excitation.lower() == 'random':
419 kwargs = {'shape': (num_signals,),
420 'n_samples': frame_length * integration_oversample * num_averages + int(steady_state_time * sample_rate),
421 'rms': excitation_level,
422 'dt': dt,
423 'low_frequency_cutoff': excitation_min_frequency,
424 'high_frequency_cutoff': bandwidth if excitation_max_frequency is None else excitation_max_frequency}
425 kwargs.update(generator_kwargs)
426 signals = generator.random(**kwargs)
427 elif excitation.lower() == 'burst random':
428 kwargs = {'shape': (num_signals,),
429 'n_samples': frame_length * integration_oversample,
430 'on_fraction': signal_fraction,
431 'delay_fraction': 0,
432 'rms': excitation_level,
433 'dt': dt,
434 'low_frequency_cutoff': excitation_min_frequency,
435 'high_frequency_cutoff': bandwidth if excitation_max_frequency is None else excitation_max_frequency}
436 kwargs.update(generator_kwargs)
437 signal_list = [generator.burst_random(**kwargs) for i in range(num_averages)]
438 full_list = []
439 for i, signal in enumerate(signal_list):
440 full_list.append(
441 np.zeros((num_signals, int(extra_time_between_frames * sample_rate))))
442 full_list.append(signal)
443 signals = np.concatenate(full_list, axis=-1)
444 elif excitation.lower() == 'chirp':
445 if num_signals > 1:
446 print('Warning: Chirp generally not recommended for multi-reference excitation.')
447 kwargs = {'frequency_min': 0 if excitation_min_frequency is None else excitation_min_frequency,
448 'frequency_max': bandwidth if excitation_max_frequency is None else excitation_max_frequency,
449 'signal_length': frame_time,
450 'dt': dt}
451 kwargs.update(generator_kwargs)
452 signals = np.array([
453 generator.chirp(**kwargs)
454 for i in range(num_signals)
455 ]) * excitation_level
456 signals = np.tile(signals, [1, num_averages +
457 int(np.ceil(steady_state_time / frame_time))])
458 elif excitation.lower() == 'hammer':
459 if num_signals > 1:
460 print(
461 'Warning: Hammer impact generally not recommended for multi-reference excitation, consider multi-hammer instead')
462 pulse_width = 2 / (bandwidth if excitation_max_frequency is None else excitation_max_frequency)
463 signal_length = int(frame_length * integration_oversample * num_averages + (num_averages + 1)
464 * extra_time_between_frames * sample_rate + 2 * pulse_width * sample_rate)
465 pulse_times = np.arange(num_averages)[
466 :, np.newaxis] * (frame_time + extra_time_between_frames) + pulse_width + extra_time_between_frames
467 kwargs = {'signal_length': signal_length,
468 'pulse_time': pulse_times,
469 'pulse_width': pulse_width,
470 'pulse_peak': excitation_level,
471 'dt': dt,
472 'sine_exponent': 2}
473 kwargs.update(generator_kwargs)
474 signals = generator.pulse(**kwargs)
475 signals = np.tile(signals, [num_signals, 1])
476 elif excitation.lower() == 'multi-hammer':
477 signal_length = frame_length * integration_oversample
478 pulse_width = 2 / (bandwidth if excitation_max_frequency is None else excitation_max_frequency)
479 signals = []
480 for i in range(num_signals):
481 signals.append([])
482 for j in range(num_averages):
483 pulse_times = []
484 last_pulse = 0
485 while last_pulse < frame_time * signal_fraction:
486 next_pulse = last_pulse + pulse_width * (np.random.rand() * multihammer_impact_spacing_factor + 1)
487 pulse_times.append(next_pulse)
488 last_pulse = next_pulse
489 pulse_times = np.array(pulse_times)
490 pulse_times = pulse_times[pulse_times <
491 frame_time * signal_fraction, np.newaxis]
492 kwargs = {'signal_length': signal_length,
493 'pulse_time': pulse_times,
494 'pulse_width': pulse_width,
495 'pulse_peak': excitation_level,
496 'dt': dt,
497 'sine_exponent': 2}
498 kwargs.update(generator_kwargs)
499 signal = generator.pulse(**kwargs)
500 signals[-1].append(np.zeros(int(extra_time_between_frames * sample_rate)))
501 signals[-1].append(signal)
502 signals[-1].append(np.zeros(int(extra_time_between_frames * sample_rate)))
503 signals[-1] = np.concatenate(signals[-1], axis=-1)
504 signals = np.array(signals)
505 elif excitation.lower() == 'sine':
506 if num_signals > 1:
507 print(
508 'Warning: Sine signal generally not recommended for multi-reference excitation')
509 frequencies = excitation_max_frequency if excitation_min_frequency is None else excitation_min_frequency
510 num_samples = frame_length * integration_oversample * num_averages + int(steady_state_time * sample_rate)
511 kwargs = {'frequencies': frequencies,
512 'dt': dt,
513 'num_samples': num_samples,
514 'amplitudes': excitation_level}
515 kwargs.update(generator_kwargs)
516 signals = np.tile(generator.sine(**kwargs), (num_signals, 1))
517 # Set up the integration
518 responses, references = self.time_integrate(
519 signals, dt, responses, references, displacement_derivative,**time_integrate_kwargs)
520 # Now add noise
521 responses.ordinate += response_noise_level * np.random.randn(*responses.ordinate.shape)
522 references.ordinate += excitation_noise_level * np.random.randn(*references.ordinate.shape)
523 # Filter with antialiasing filters, divide filter order by 2 because of filtfilt
524 if antialias_filter_order > 0:
525 lowpass_b, lowpass_a = butter(antialias_filter_order // 2,
526 antialias_filter_cutoff_factor * bandwidth, fs=sample_rate)
527 responses.ordinate = filtfilt(lowpass_b, lowpass_a, responses.ordinate)
528 references.ordinate = filtfilt(lowpass_b, lowpass_a, references.ordinate)
529 if integration_oversample > 1:
530 responses = responses.downsample(integration_oversample)
531 references = references.downsample(integration_oversample)
532 responses = responses.extract_elements_by_abscissa(steady_state_time, np.inf)
533 references = references.extract_elements_by_abscissa(steady_state_time, np.inf)
534 return responses, references
536 def copy(self):
537 return SystemNL(self.coordinate.copy(),self.state_derivative_function,
538 self.state_size,self.displacement_output_function,self.velocity_output_function,
539 self.acceleration_output_function,self.excitation_output_function)