From ca9d4f4304fb733d4096311c409c6585dcbebf2a Mon Sep 17 00:00:00 2001 From: Ilhan Polat Date: Fri, 21 Sep 2018 11:35:14 +0200 Subject: [PATCH 1/2] FIX: step_resp : Added missing integrator counts Array sizes didn't match since integrators were not taken into account. Discovered via beam.mat from SLICOT collections --- harold/_time_domain.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/harold/_time_domain.py b/harold/_time_domain.py index 96fcff2..f4b8e61 100644 --- a/harold/_time_domain.py +++ b/harold/_time_domain.py @@ -412,8 +412,8 @@ def _compute_tfinal_and_dt(sys, is_step=True): # The rest texp_mode = np.log(decay_percent) / np.abs(psub[~iw & ~ints].real) tfinal += texp_mode.tolist() - dt += minimum(texp_mode / 50, (2 * np.pi / pts_per_cycle / wnsub[~iw]) - ).tolist() + dt += minimum(texp_mode / 50, + (2 * np.pi / pts_per_cycle / wnsub[~iw & ~ints])).tolist() # All integrators? if len(tfinal) == 0: From 99dfc7e0c3fec34fda5c2f86c63094ebbaffe5b6 Mon Sep 17 00:00:00 2001 From: Ilhan Polat Date: Wed, 10 Oct 2018 08:05:33 +0200 Subject: [PATCH 2/2] MAINT: Temporary fix for continuous case --- harold/_time_domain.py | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/harold/_time_domain.py b/harold/_time_domain.py index f4b8e61..29e2eda 100644 --- a/harold/_time_domain.py +++ b/harold/_time_domain.py @@ -324,7 +324,7 @@ def _compute_tfinal_and_dt(sys, is_step=True): default_tfinal = 5 # Default simulation horizon total_cycles = 5 # number of cycles for oscillating modes pts_per_cycle = 25 # Number of points divide a period of oscillation - decay_percent = 100 # Factor of reduction for real pole decays + log_decay_percent = np.log(100) # Factor of reduction for real pole decays # if a static model is given, don't bother with checks if sys._isgain: @@ -335,16 +335,24 @@ def _compute_tfinal_and_dt(sys, is_step=True): if sys._isdiscrete: # System already has sampling fixed hence we can't fall into the same - # trap mentioned above. Just get nonintegrating slow modes. + # trap mentioned above. Just get nonintegrating slow modes together + # with the damping. dt = sys._dt p = eigvals(sys.a) ps = np.log(p)/sys._dt ps[np.abs(ps) <= sqrt_eps] = 0. wn = np.abs(ps) + zeta = -np.cos(np.angle(ps)) - if np.any(ps.real < 0): - stable_decay = np.max(ps.real < 0.) - tfinal = np.log(decay_percent) / np.abs(stable_decay) + if np.any(ps.real < sqrt_eps): + # Let the extra underdamped ones ring + ind = (ps.real < 0) & (zeta > 1e-3) + + stable_decay = np.abs(wn[ind] * zeta[ind]).max() + stable_decay = np.reciprocal(stable_decay) + tfinal = np.min([max_points_z * sys._dt, + log_decay_percent / stable_decay]) + tfinal = np.ceil(tfinal / sys._dt) * sys._dt else: tfinal = dt * min_points_z @@ -358,7 +366,6 @@ def _compute_tfinal_and_dt(sys, is_step=True): elif np.any(ps.real > 0): tfinal = tfinal*2 - dt = tfinal / max_points_z if tfinal // dt > max_points_z else dt tfinal = dt*min_points_z if tfinal // dt < min_points_z else tfinal return tfinal, dt @@ -409,8 +416,8 @@ def _compute_tfinal_and_dt(sys, is_step=True): if np.any(iw): tfinal += (total_cycles * 2 * np.pi / wnsub[iw]).tolist() dt += (2 * np.pi / pts_per_cycle / wnsub[iw]).tolist() - # The rest - texp_mode = np.log(decay_percent) / np.abs(psub[~iw & ~ints].real) + # The rest ~ts = log(%ss value) / exp(Re(λ)t) + texp_mode = log_decay_percent / np.abs(psub[~iw & ~ints].real) tfinal += texp_mode.tolist() dt += minimum(texp_mode / 50, (2 * np.pi / pts_per_cycle / wnsub[~iw & ~ints])).tolist()