allow solve_ivp errors and improve arrow placement

murrayrm · murrayrm · commit d44a5771e4b7 · 2024-06-04T22:15:27.000-07:00
diff --git a/control/freqplot.py b/control/freqplot.py
@@ -1963,9 +1963,22 @@ def _add_arrows_to_line2D(
     if transform is None:
         transform = axes.transData
 
+    # Figure out the size of the axes (length of diagonal)
+    xlim, ylim = axes.get_xlim(), axes.get_ylim()
+    ul, lr = np.array([xlim[0], ylim[0]]), np.array([xlim[1], ylim[1]])
+    diag = np.linalg.norm(ul - lr)
+
     # Compute the arc length along the curve
     s = np.cumsum(np.sqrt(np.diff(x) ** 2 + np.diff(y) ** 2))
 
+    # Truncate the number of arrows if the curve is short
+    # TODO: figure out a smarter way to do this
+    frac = min(s[-1] / diag, 1)
+    if len(arrow_locs) and frac < 0.05:
+        arrow_locs = []         # too short; no arrows at all
+    elif len(arrow_locs) and frac < 0.2:
+        arrow_locs = [0.5]      # single arrow in the middle
+
     arrows = []
     for loc in arrow_locs:
         n = np.searchsorted(s, s[-1] * loc)
diff --git a/control/nlsys.py b/control/nlsys.py
@@ -1317,7 +1317,7 @@ def nlsys(
 
 
 def input_output_response(
-        sys, T, U=0., X0=0, params=None,
+        sys, T, U=0., X0=0, params=None, ignore_error=False,
         transpose=False, return_x=False, squeeze=None,
         solve_ivp_kwargs=None, t_eval='T', **kwargs):
     """Compute the output response of a system to a given input.
@@ -1593,7 +1593,7 @@ def ivp_rhs(t, x):
         soln = sp.integrate.solve_ivp(
             ivp_rhs, (T0, Tf), X0, t_eval=t_eval,
             vectorized=False, **solve_ivp_kwargs)
-        if not soln.success:
+        if not ignore_error and not soln.success:
             raise RuntimeError("solve_ivp failed: " + soln.message)
 
         # Compute inputs and outputs for each time point
diff --git a/control/phaseplot.py b/control/phaseplot.py
@@ -375,7 +375,7 @@ def streamlines(
             sys, revsys, timepts, X0, params, dir,
             gridtype=gridtype, gridspec=gridspec, xlim=xlim, ylim=ylim)
 
-        # Plot the trajectory
+        # Plot the trajectory (if there is one)
         if traj.shape[1] > 1:
             out.append(
                 ax.plot(traj[0], traj[1], color=color))
@@ -596,6 +596,7 @@ def separatrices(
                     color = unstable_color
                     linestyle = '-'
 
+                # Plot the trajectory (if there is one)
                 if traj.shape[1] > 1:
                     out.append(ax.plot(
                         traj[0], traj[1], color=color, linestyle=linestyle))
@@ -883,12 +884,13 @@ def _create_trajectory(
         gridtype=None, gridspec=None, xlim=None, ylim=None):
     # Comput ethe forward trajectory
     if dir == 'forward' or dir == 'both':
-        fwdresp = input_output_response(sys, timepts, X0=X0, params=params)
+        fwdresp = input_output_response(
+            sys, timepts, X0=X0, params=params, ignore_error=True)
 
     # Compute the reverse trajectory
     if dir == 'reverse' or dir == 'both':
         revresp = input_output_response(
-            revsys, timepts, X0=X0, params=params)
+            revsys, timepts, X0=X0, params=params, ignore_error=True)
 
     # Create the trace to plot
     if dir == 'forward':
@@ -898,7 +900,14 @@ def _create_trajectory(
     elif dir == 'both':
         traj = np.hstack([revresp.states[:, :1:-1], fwdresp.states])
 
-    return traj
+    # Remove points outside the window (keep first point beyond boundary)
+    inrange = np.asarray(
+        (traj[0] >= xlim[0]) & (traj[0] <= xlim[1]) &
+        (traj[1] >= ylim[0]) & (traj[1] <= ylim[1]))
+    inrange[:-1] = inrange[:-1] | inrange[1:]   # keep if next point in range
+    inrange[1:] = inrange[1:] | inrange[:-1]    # keep if prev point in range
+
+    return traj[:, inrange]
 
 
 def _make_timepts(timepts, i):
