Merge pull request #1051 from pymor/bode_plot

Add Bode plot

docs/source/tutorial_bt.rst

@@ -54,7 +54,7 @@ In particular, we will use the
 :meth:`~pymor.models.iosys.LTIModel.from_matrices` method of |LTIModel|, which
 instantiates an |LTIModel| from NumPy or SciPy matrices.
 
-First, we do the necessary imports.
+First, we do the necessary imports and some matplotlib style choices.
 
 .. jupyter-execute::
 
@@ -64,6 +64,8 @@ First, we do the necessary imports.
     from pymor.models.iosys import LTIModel
     from pymor.reductors.bt import BTReductor
 
+    plt.rcParams['axes.grid'] = True
+
 Next, we can assemble the matrices based on a centered finite difference
 approximation:
 
@@ -112,14 +114,23 @@ We can also see some basic information from `fom`'s string representation
     print(fom)
 
 To visualize the behavior of the `fom`, we can draw its magnitude plot, i.e.,
-a visualization of the mapping :math:`\omega \mapsto H(i \omega)`, where
+a visualization of the mapping :math:`\omega \mapsto H(\imath \omega)`, where
 :math:`H(s) = C (s E - A)^{-1} B + D` is the transfer function of the system.
 
 .. jupyter-execute::
 
     w = np.logspace(-2, 8, 50)
-    fom.mag_plot(w)
-    plt.grid()
+    _ = fom.mag_plot(w)
+
+We can also see the Bode plot, which shows the magnitude and phase of the
+components of the transfer function.
+In particular, :math:`\lvert H_{ij}(\imath \omega) \rvert` is in subplot
+:math:`(2 i - 1, j)` and :math:`\arg(H_{ij}(\imath \omega))` is in subplot
+:math:`(2 i, j)`.
+
+.. jupyter-execute::
+
+    _ = fom.bode_plot(w)
 
 Plotting the Hankel singular values shows us how well an LTI system can be
 approximated by a reduced-order model
@@ -129,8 +140,7 @@ approximated by a reduced-order model
     hsv = fom.hsv()
     fig, ax = plt.subplots()
     ax.semilogy(range(1, len(hsv) + 1), hsv, '.-')
-    ax.set_title('Hankel singular values')
-    ax.grid()
+    _ = ax.set_title('Hankel singular values')
 
 As expected for a heat equation, the Hankel singular values decay rapidly.
 
@@ -218,8 +228,7 @@ based on the Hankel singular values:
     ax.semilogy(range(1, len(error_bounds) + 1), error_bounds, '.-')
     ax.semilogy(range(1, len(hsv)), hsv[1:], '.-')
     ax.set_xlabel('Reduced order')
-    ax.set_title(r'Upper and lower $\mathcal{H}_\infty$ error bounds')
-    ax.grid()
+    _ = ax.set_title(r'Upper and lower $\mathcal{H}_\infty$ error bounds')
 
 To get a reduced-order model of order 10, we call the `reduce` method with the
 appropriate argument:
@@ -241,23 +250,36 @@ models
     fig, ax = plt.subplots()
     fom.mag_plot(w, ax=ax, label='FOM')
     rom.mag_plot(w, ax=ax, linestyle='--', label='ROM')
-    ax.legend()
-    ax.grid()
+    _ = ax.legend()
+
+as well as Bode plots
+
+.. jupyter-execute::
+
+    fig, axs = plt.subplots(6, 2, figsize=(12, 24), sharex=True, constrained_layout=True)
+    fom.bode_plot(w, ax=axs)
+    _ = rom.bode_plot(w, ax=axs, linestyle='--')
+
+Also, we can plot the magnitude plot of the error system
+
+.. jupyter-execute::
+
+    err = fom - rom
+    _ = err.mag_plot(w)
 
-and plot the magnitude plot of the error system
+and its Bode plot
 
 .. jupyter-execute::
 
-    (fom - rom).mag_plot(w)
-    plt.grid()
+    _ = err.bode_plot(w)
 
 We can compute the relative errors in :math:`\mathcal{H}_\infty` or
 :math:`\mathcal{H}_2` (or Hankel) norm
 
 .. jupyter-execute::
 
-    print(f'Relative Hinf error: {(fom - rom).hinf_norm() / fom.hinf_norm():.3e}')
-    print(f'Relative H2 error:   {(fom - rom).h2_norm() / fom.h2_norm():.3e}')
+    print(f'Relative Hinf error: {err.hinf_norm() / fom.hinf_norm():.3e}')
+    print(f'Relative H2 error:   {err.h2_norm() / fom.h2_norm():.3e}')
 
 .. note::
 

src/pymor/models/iosys.py

@@ -16,8 +16,7 @@
                                    SecondOrderModelOperator)
 from pymor.operators.constructions import IdentityOperator, LincombOperator, ZeroOperator
 from pymor.operators.numpy import NumpyMatrixOperator
-from pymor.parameters.base import Mu, Parameters
-from pymor.tools.formatrepr import indent_value
+from pymor.parameters.base import Mu
 from pymor.vectorarrays.block import BlockVectorSpace
 
 
@@ -81,6 +80,97 @@ def freq_resp(self, w, mu=None):
         assert self.parameters.assert_compatible(mu)
         return np.stack([self.eval_tf(1j * wi, mu=mu) for wi in w])
 
+    def bode(self, w, mu=None):
+        """Compute magnitudes and phases.
+
+        Parameters
+        ----------
+        w
+            A sequence of angular frequencies at which to compute the transfer function.
+        mu
+            |Parameter values| for which to evaluate the transfer function.
+
+        Returns
+        -------
+        mag
+            Transfer function magnitudes at frequencies in `w`, |NumPy array| of shape
+            `(len(w), self.output_dim, self.input_dim)`.
+        phase
+            Transfer function phases (in radians) at frequencies in `w`, |NumPy array| of shape
+            `(len(w), self.output_dim, self.input_dim)`.
+        """
+        w = np.asarray(w)
+        mag = np.abs(self.freq_resp(w, mu=mu))
+        phase = np.angle(self.freq_resp(w, mu=mu))
+        phase = np.unwrap(phase, axis=0)
+        return mag, phase
+
+    def bode_plot(self, w, mu=None, ax=None, Hz=False, dB=False, deg=True, **mpl_kwargs):
+        """Draw the Bode plot for all input-output pairs.
+
+        Parameters
+        ----------
+        w
+            A sequence of angular frequencies at which to compute the transfer function.
+        mu
+            |Parameter| for which to evaluate the transfer function.
+        ax
+            Axis of shape (2 * `self.output_dim`, `self.input_dim`) to which to plot.
+            If not given, `matplotlib.pyplot.gcf` is used to get the figure and create axis.
+        Hz
+            Should the frequency be in Hz on the plot.
+        dB
+            Should the magnitude be in dB on the plot.
+        deg
+            Should the phase be in degrees (otherwise in radians).
+        mpl_kwargs
+            Keyword arguments used in the matplotlib plot function.
+
+        Returns
+        -------
+        artists
+            List of matplotlib artists added.
+        """
+        if ax is None:
+            import matplotlib.pyplot as plt
+            fig = plt.gcf()
+            width, height = plt.rcParams['figure.figsize']
+            fig.set_size_inches(self.input_dim * width, 2 * self.output_dim * height)
+            fig.set_constrained_layout(True)
+            ax = fig.subplots(2 * self.output_dim, self.input_dim, sharex=True, squeeze=False)
+        else:
+            assert isinstance(ax, np.ndarray) and ax.shape == (2 * self.output_dim, self.input_dim)
+            fig = ax[0, 0].get_figure()
+
+        w = np.asarray(w)
+        freq = w / (2 * np.pi) if Hz else w
+        mag, phase = self.bode(w, mu=mu)
+        if deg:
+            phase *= 180 / np.pi
+
+        artists = np.empty_like(ax)
+        freq_label = f'Frequency ({"Hz" if Hz else "rad/s"})'
+        mag_label = f'Magnitude{" (dB)" if dB else ""}'
+        phase_label = f'Phase ({"deg" if deg else "rad"})'
+        for i in range(self.output_dim):
+            for j in range(self.input_dim):
+                if dB:
+                    artists[2 * i, j] = ax[2 * i, j].semilogx(freq, 20 * np.log2(mag[:, i, j]),
+                                                              **mpl_kwargs)
+                else:
+                    artists[2 * i, j] = ax[2 * i, j].loglog(freq, mag[:, i, j],
+                                                            **mpl_kwargs)
+                artists[2 * i + 1, j] = ax[2 * i + 1, j].semilogx(freq, phase[:, i, j],
+                                                                  **mpl_kwargs)
+        for i in range(self.output_dim):
+            ax[2 * i, 0].set_ylabel(mag_label)
+            ax[2 * i + 1, 0].set_ylabel(phase_label)
+        for j in range(self.input_dim):
+            ax[-1, j].set_xlabel(freq_label)
+        fig.suptitle('Bode plot')
+
+        return artists
+
     def mag_plot(self, w, mu=None, ax=None, ord=None, Hz=False, dB=False, **mpl_kwargs):
         """Draw the magnitude plot.