.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/ode_dynamics/stiffness/plot_01_dahlquist_a_stability.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_ode_dynamics_stiffness_plot_01_dahlquist_a_stability.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_ode_dynamics_stiffness_plot_01_dahlquist_a_stability.py:


Dahlquist's A-stability: backward Euler vs RK4
==============================================

Applied to the test equation :math:`y' = \lambda y`, a one-step method
multiplies :math:`y` by a stability function :math:`R(z)`, with
:math:`z = \lambda h`. Dahlquist (1963) called a method *A-stable* if
:math:`|R(z)| \le 1` on the whole left half-plane, so that no decaying
mode can ever be amplified, whatever the step size. Explicit methods
never qualify, since their :math:`R` is a polynomial and grows without
bound. RK4's region reaches only :math:`z \approx -2.785` on the real
axis. Backward Euler's :math:`R(z) = 1/(1 - z)` is A-stable. Dahlquist
also proved that no A-stable linear multistep method can exceed order 2
(his "second barrier").

.. GENERATED FROM PYTHON SOURCE LINES 18-24

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from mathematicskit.integrators import implicit_euler_integrate, is_absolutely_stable, njit, rk4_integrate
    from mathematicskit.ode_dynamics.visualizers import plot_stability_regions








.. GENERATED FROM PYTHON SOURCE LINES 25-27

Stability regions in the complex z-plane
----------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 27-39

.. code-block:: Python


    fig, axes = plt.subplots(1, 3, figsize=(14, 4.2))
    for ax, (method, name) in zip(axes[:2], [("rk4", "RK4"), ("implicit_euler", "backward Euler")], strict=True):
        plot_stability_regions(method, ax=ax, re_range=(-5, 3), im_range=(-4, 4), labels=[name])
        ax.set_title(f"{name}: shaded where |R(z)| <= 1")

    x, y = np.meshgrid(np.linspace(-50, 0, 201), np.linspace(-50, 50, 401))
    left_half_plane = x + 1j * y
    for method in ("rk4", "implicit_euler"):
        print(f"{method}: stable on the whole sampled left half-plane? {is_absolutely_stable(method, left_half_plane).all()}")





.. image-sg:: /api/gallery/ode_dynamics/stiffness/images/sphx_glr_plot_01_dahlquist_a_stability_001.png
   :alt: RK4: shaded where |R(z)| <= 1, backward Euler: shaded where |R(z)| <= 1
   :srcset: /api/gallery/ode_dynamics/stiffness/images/sphx_glr_plot_01_dahlquist_a_stability_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none

    rk4: stable on the whole sampled left half-plane? False
    implicit_euler: stable on the whole sampled left half-plane? True




.. GENERATED FROM PYTHON SOURCE LINES 40-46

Consequence on a stiff equation, lambda * h = 30
------------------------------------------------

The Prothero-Robinson equation :math:`y' = -\lambda(y - \cos t) - \sin t`
has the smooth solution :math:`y = \cos t`, but with :math:`\lambda = 3000`
and :math:`h = 0.01`, :math:`z = -30` lies far outside RK4's region.

.. GENERATED FROM PYTHON SOURCE LINES 46-67

.. code-block:: Python



    @njit
    def prothero_robinson(state, t, params):
        return -params[0] * (state - np.cos(t)) - np.sin(t)


    params = np.array([3000.0])
    ts, ys_ie = implicit_euler_integrate(prothero_robinson, np.array([1.0]), 0.0, 1e-2, 40, params)
    _, ys_rk4 = rk4_integrate(prothero_robinson, np.array([1.0]), 0.0, 1e-2, 40, params)
    axes[2].semilogy(ts, np.abs(ys_rk4[:, 0] - np.cos(ts)) + 1e-16, "o-", ms=3, label="RK4")
    axes[2].semilogy(ts, np.abs(ys_ie[:, 0] - np.cos(ts)) + 1e-16, "s-", ms=3, label="backward Euler")
    axes[2].set_xlabel("t")
    axes[2].set_ylabel("|error|")
    axes[2].set_title("z = -30: only the A-stable method survives")
    axes[2].legend()
    fig.tight_layout()

    print(f"final error: RK4 {abs(ys_rk4[-1, 0] - np.cos(ts[-1])):.2e}, backward Euler {abs(ys_ie[-1, 0] - np.cos(ts[-1])):.2e}")

    plt.show()




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

 .. code-block:: none

    final error: RK4 7.41e+172, backward Euler 1.54e-06





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

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


.. _sphx_glr_download_api_gallery_ode_dynamics_stiffness_plot_01_dahlquist_a_stability.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/ode_dynamics/stiffness/plot_01_dahlquist_a_stability.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

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

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

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

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

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

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


.. only:: html

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

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