.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/signal_processing/fourier/plot_01_gibbs_phenomenon.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` 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_signal_processing_fourier_plot_01_gibbs_phenomenon.py: Gibbs's phenomenon (1899) ========================= Albert Michelson's harmonic analyser (1898) summed 80 Fourier terms mechanically, and its curves for a square wave overshot at every jump. Michelson suspected the machine. Josiah Willard Gibbs (1899) showed in *Nature* that the mathematics itself overshoots: near a discontinuity the partial sums (:func:`~engineeringkit.signal_processing.fourier_partial_sum`) rise about 9% of the jump above it, however many terms are kept. The overshoot only narrows toward the jump. Henry Wilbraham had found the same effect in 1848, unnoticed. .. GENERATED FROM PYTHON SOURCE LINES 16-22 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from scipy.special import sici from engineeringkit.signal_processing import fourier_partial_sum .. GENERATED FROM PYTHON SOURCE LINES 23-28 A square wave and its partial sums ---------------------------------- The square wave of period 1 and amplitude :math:`\pm 1` has :math:`b_k = 4/(\pi k)` for odd :math:`k` and no other terms. .. GENERATED FROM PYTHON SOURCE LINES 28-49 .. code-block:: Python period = 1.0 t = np.linspace(-0.25, 0.25, 20001) limit = 2 / np.pi * sici(np.pi)[0] # peak of the partial sums as K grows print(f"predicted limit of the peak: {limit:.5f} (overshoot {(limit - 1) / 2:.2%} of the jump)") fig, axes = plt.subplots(1, 2, figsize=(10, 3.8)) peaks = {} for K in (5, 25, 125, 625): k = np.arange(1, K + 1) b = np.where(k % 2 == 1, 4 / (np.pi * k), 0.0) x = fourier_partial_sum(t, period, np.zeros(K), b) peaks[K] = x.max() position = t[np.argmax(x)] print(f"K = {K:3d}: peak {peaks[K]:.5f} at t = {position:.5f} s") if K <= 125: axes[0].plot(t, x, lw=0.8, label=f"K = {K}") axes[0].plot(t, np.sign(t), "k", lw=1.5, label="square wave") axes[0].set(xlabel="t [s]", ylabel="x", title="Partial sums near the jump at t = 0") axes[0].legend(fontsize=8) .. image-sg:: /api/gallery/signal_processing/fourier/images/sphx_glr_plot_01_gibbs_phenomenon_001.png :alt: Partial sums near the jump at t = 0 :srcset: /api/gallery/signal_processing/fourier/images/sphx_glr_plot_01_gibbs_phenomenon_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none predicted limit of the peak: 1.17898 (overshoot 8.95% of the jump) K = 5: peak 1.18836 at t = 0.08333 s K = 25: peak 1.17947 at t = 0.01922 s K = 125: peak 1.17900 at t = 0.00398 s K = 625: peak 1.17898 at t = 0.00080 s .. GENERATED FROM PYTHON SOURCE LINES 50-55 The overshoot does not shrink ----------------------------- The peak tends to :math:`\frac{2}{\pi}\operatorname{Si}(\pi) \approx 1.179`, not to 1. It sits at :math:`t \approx T/(2K)`, moving toward the jump. .. GENERATED FROM PYTHON SOURCE LINES 55-59 .. code-block:: Python assert abs(peaks[625] - limit) < 1e-3 assert all(p > 1.17 for p in peaks.values()) .. GENERATED FROM PYTHON SOURCE LINES 60-65 What does converge: the energy of the error ------------------------------------------- The mean squared error over a period falls like :math:`1/K`. Fourier series converge in energy even where they never converge uniformly. .. GENERATED FROM PYTHON SOURCE LINES 65-83 .. code-block:: Python t_period = np.linspace(-0.5, 0.5, 40001)[:-1] Ks = np.array([5, 25, 125, 625]) mse = [] for K in Ks: k = np.arange(1, K + 1) b = np.where(k % 2 == 1, 4 / (np.pi * k), 0.0) mse.append(np.mean((fourier_partial_sum(t_period, period, np.zeros(K), b) - np.sign(t_period)) ** 2)) slope = np.polyfit(np.log(Ks), np.log(mse), 1)[0] print(f"mean squared error falls as K^{slope:.2f}") assert abs(slope + 1) < 0.1 axes[1].loglog(Ks, mse, "o-", label="mean squared error") axes[1].loglog(Ks, np.array([peaks[K] for K in Ks]) - 1, "s-", label="overshoot") axes[1].set(xlabel="number of terms K", title="Error energy vanishes, overshoot does not") axes[1].legend() fig.tight_layout() plt.show() .. rst-class:: sphx-glr-script-out .. code-block:: none mean squared error falls as K^-0.97 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.431 seconds) .. _sphx_glr_download_api_gallery_signal_processing_fourier_plot_01_gibbs_phenomenon.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/signal_processing/fourier/plot_01_gibbs_phenomenon.ipynb :alt: Launch JupyterLite :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_01_gibbs_phenomenon.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_01_gibbs_phenomenon.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_01_gibbs_phenomenon.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_