.. 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_04_radau_iia.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_04_radau_iia.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_04_radau_iia.py:


Butcher, Ehle and Radau IIA: L-stable implicit Runge-Kutta
==========================================================

Butcher (1964) built implicit Runge-Kutta methods on Gauss and Radau
quadrature nodes. Ehle (1969) showed that their stability functions are
Padé approximants of :math:`e^z`, and that the Radau IIA family is
A-stable and, moreover, *L-stable*: :math:`R(z) \to 0` as
:math:`z \to -\infty`, so infinitely stiff modes are damped in a single
step. The three-stage, order-5 method used by
:func:`~mathematicskit.integrators.stiff_integrate` has

.. math::

    R(z) = \frac{1 + \tfrac25 z + \tfrac1{20} z^2}
                {1 - \tfrac35 z + \tfrac3{20} z^2 - \tfrac1{60} z^3}.

The trapezoidal rule is A-stable but not L-stable, since
:math:`R(z) \to -1`: stiff components flip sign each step instead of dying.

.. GENERATED FROM PYTHON SOURCE LINES 23-28

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from mathematicskit.integrators import njit, stiff_integrate








.. GENERATED FROM PYTHON SOURCE LINES 29-31

Damping of stiff modes
----------------------

.. GENERATED FROM PYTHON SOURCE LINES 31-46

.. code-block:: Python


    z = -np.logspace(-1, 4, 300)
    r_radau = (1 + 2 * z / 5 + z**2 / 20) / (1 - 3 * z / 5 + 3 * z**2 / 20 - z**3 / 60)
    r_trap = (1 + z / 2) / (1 - z / 2)

    fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(11, 4.2))
    ax0.loglog(-z, np.abs(r_radau), label="Radau IIA (order 5)")
    ax0.loglog(-z, np.abs(r_trap), label="trapezoidal rule (order 2)")
    ax0.loglog(-z, np.exp(z).clip(1e-30), "k:", label=r"exact $e^{z}$")
    ax0.set_ylim(1e-12, 2)
    ax0.set_xlabel(r"$-z = |\lambda| h$")
    ax0.set_ylabel(r"$|R(z)|$")
    ax0.set_title("L-stability: R(z) -> 0 as z -> -infinity")
    ax0.legend()




.. image-sg:: /api/gallery/ode_dynamics/stiffness/images/sphx_glr_plot_04_radau_iia_001.png
   :alt: L-stability: R(z) -> 0 as z -> -infinity
   :srcset: /api/gallery/ode_dynamics/stiffness/images/sphx_glr_plot_04_radau_iia_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none


    <matplotlib.legend.Legend object at 0x7fd7c3149d00>



.. GENERATED FROM PYTHON SOURCE LINES 47-54

The stiff van der Pol oscillator, mu = 1000
-------------------------------------------

Relaxation oscillations with slow drifts and near-instant jumps. This is
the classic stiff Radau IIA benchmark (Hairer and Wanner). Over
:math:`0 \le t \le 3000` Radau IIA needs about 1200 steps; scipy's
explicit ``RK45`` needs about 1.7 million at the same tolerance.

.. GENERATED FROM PYTHON SOURCE LINES 54-73

.. code-block:: Python



    @njit
    def van_der_pol(state, t, params):
        mu = params[0]
        return np.array([state[1], mu * (1 - state[0] ** 2) * state[1] - state[0]])


    ts, ys = stiff_integrate(van_der_pol, np.array([2.0, 0.0]), 0.0, 3000.0, np.array([1000.0]), method="Radau", rtol=1e-6, atol=1e-8)
    ax1.plot(ts, ys[:, 0], "-", lw=1)
    ax1.plot(ts, ys[:, 0], ".", ms=2)
    ax1.set_xlabel("t")
    ax1.set_ylabel("x")
    ax1.set_title(f"Radau IIA on van der Pol, mu = 1000 ({len(ts) - 1} steps)")
    fig.tight_layout()

    print(f"Radau IIA: {len(ts) - 1} steps; smallest step {np.min(np.diff(ts)):.1e}, largest {np.max(np.diff(ts)):.1e}")

    plt.show()




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

 .. code-block:: none

    Radau IIA: 1187 steps; smallest step 3.9e-05, largest 1.1e+02





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

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


.. _sphx_glr_download_api_gallery_ode_dynamics_stiffness_plot_04_radau_iia.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_04_radau_iia.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

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

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

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

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

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

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


.. only:: html

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

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