.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/mechanical_engineering/vibration/plot_04_rayleigh_forced_vibration.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_api_gallery_mechanical_engineering_vibration_plot_04_rayleigh_forced_vibration.py>`
        to download the full example code or to run this example in your browser via JupyterLite.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_api_gallery_mechanical_engineering_vibration_plot_04_rayleigh_forced_vibration.py:


Rayleigh's *Theory of Sound*: forced vibration and resonance (1877)
===================================================================

The third chapter of Rayleigh's *Theory of Sound* drives the damped
oscillator with a harmonic force, :math:`m\ddot x + c\dot x + kx =
F_0\cos\omega t`. Once the free motion has died away the mass moves at the
forcing frequency with amplitude :math:`(F_0/k)\,M`, the static deflection
times the magnification factor
(:func:`~engineeringkit.mechanical_engineering.magnification_factor`), and
lags the force by a phase
(:func:`~engineeringkit.mechanical_engineering.phase_angle`) that passes
through 90 degrees at :math:`\omega = \omega_n` whatever the damping. Below
resonance the spring carries the force, at resonance only damping limits
the motion, and above it the mass's inertia does. The example compares a
simulated run-up of three forcing frequencies with the closed form.

.. GENERATED FROM PYTHON SOURCE LINES 20-25

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from engineeringkit.mechanical_engineering import magnification_factor, phase_angle, simulate_sdof








.. GENERATED FROM PYTHON SOURCE LINES 26-31

Steady state of a simulated oscillator
--------------------------------------

After the transient decays, fit the amplitude and phase of the simulated
displacement and compare them with :math:`M` and :math:`\phi`.

.. GENERATED FROM PYTHON SOURCE LINES 31-49

.. code-block:: Python


    m, k, c, F0 = 4.0, 1.6e4, 40.0, 100.0
    wn, zeta = np.sqrt(k / m), c / (2 * np.sqrt(k * m))
    t = np.linspace(0, 6, 60001)
    tail = t > 4
    for r in (0.5, 1.0, 2.0):
        w = r * wn
        x = simulate_sdof(m, c, k, t=t, force=lambda s, w=w: F0 * np.cos(w * s)).displacement
        basis = np.column_stack((np.cos(w * t[tail]), np.sin(w * t[tail])))
        a, b = np.linalg.lstsq(basis, x[tail], rcond=None)[0]
        amplitude, lag = np.hypot(a, b) / (F0 / k), np.arctan2(b, a)
        print(
            f"r = {r}: M = {amplitude:.4f} (closed form {magnification_factor(r, zeta):.4f}), "
            f"lag = {np.degrees(lag):.2f} deg (closed form {np.degrees(phase_angle(r, zeta)):.2f})"
        )
        assert np.isclose(amplitude, magnification_factor(r, zeta), rtol=1e-4)
        assert np.isclose(lag, phase_angle(r, zeta), atol=1e-4)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    r = 0.5: M = 1.3260 (closed form 1.3260), lag = 6.02 deg (closed form 6.02)
    r = 1.0: M = 6.3246 (closed form 6.3246), lag = 90.00 deg (closed form 90.00)
    r = 2.0: M = 0.3315 (closed form 0.3315), lag = 173.98 deg (closed form 173.98)




.. GENERATED FROM PYTHON SOURCE LINES 50-55

Resonance curves
----------------

The peak, :math:`1/(2\zeta\sqrt{1 - \zeta^2})`, sits slightly below
:math:`r = 1`; every phase curve crosses 90 degrees at :math:`r = 1`.

.. GENERATED FROM PYTHON SOURCE LINES 55-70

.. code-block:: Python


    r = np.linspace(0, 3, 3001)
    fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True, figsize=(6, 6))
    for z in (0.05, 0.1, 0.25, 0.5, 1.0):
        ax1.plot(r, magnification_factor(r, z), label=f"zeta = {z}")
        ax2.plot(r, np.degrees(phase_angle(r, z)))
        if z < 1 / np.sqrt(2):
            assert np.isclose(magnification_factor(r, z).max(), 1 / (2 * z * np.sqrt(1 - z * z)), rtol=1e-3)
        assert np.isclose(phase_angle(1.0, z), np.pi / 2)
    ax1.set(ylabel="magnification factor", ylim=(0, 11))
    ax1.legend(fontsize=8)
    ax2.set(xlabel="forcing / natural frequency", ylabel="phase lag [deg]")
    fig.tight_layout()

    plt.show()



.. image-sg:: /api/gallery/mechanical_engineering/vibration/images/sphx_glr_plot_04_rayleigh_forced_vibration_001.png
   :alt: plot 04 rayleigh forced vibration
   :srcset: /api/gallery/mechanical_engineering/vibration/images/sphx_glr_plot_04_rayleigh_forced_vibration_001.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.998 seconds)


.. _sphx_glr_download_api_gallery_mechanical_engineering_vibration_plot_04_rayleigh_forced_vibration.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: lite-badge

      .. image:: images/jupyterlite_badge_logo.svg
        :target: ../../../../lite/lab/index.html?path=api/gallery/mechanical_engineering/vibration/plot_04_rayleigh_forced_vibration.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_04_rayleigh_forced_vibration.ipynb <plot_04_rayleigh_forced_vibration.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_04_rayleigh_forced_vibration.py <plot_04_rayleigh_forced_vibration.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_04_rayleigh_forced_vibration.zip <plot_04_rayleigh_forced_vibration.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_