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

.. only:: html

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

        :ref:`Go to the end <sphx_glr_download_examples_ode_solver.py>`
        to download the full example code.

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

.. _sphx_glr_examples_ode_solver.py:


Floating-point ODE solver
=========================

This illustrates the effect of solving an ODE with forward Euler
using different floating-point formats and different quantization schemes.

The ODE in question is

.. math::

   u'(t) = v(t), \quad u(0) = 1,

   v'(t) = -u(t), \quad v(0) = 0,

and the solution is the unit circle in the (u, v) plane.

(Example was inspired from M. Croci, M. Fasi, N. J Higham, T. Mary, M. Mikaitis. "Stochastic Rounding: Implementation, Error Analysis, and Applications", in *Royal Society Open Science*, 2022.)

.. GENERATED FROM PYTHON SOURCE LINES 20-66



.. image-sg:: /examples/images/sphx_glr_ode_solver_001.png
   :alt: ode solver
   :srcset: /examples/images/sphx_glr_ode_solver_001.png
   :class: sphx-glr-single-img





.. code-block:: Python


    import math

    import matplotlib.pyplot as plt

    import apytypes as apy


    def solve_ode(exp_bits, man_bits, bias, iterations=2**14):
        h = 2 * math.pi / iterations
        ode = apy.fp([[0, h], [-h, 0]], exp_bits, man_bits, bias)
        current_value = apy.fp([1, 0], exp_bits, man_bits, bias)
        results = apy.zeros(
            (iterations + 1, 2), exp_bits=exp_bits, man_bits=man_bits, bias=bias
        )

        results[0, :] = current_value
        for i in range(iterations):
            current_value += ode @ current_value
            results[i + 1, :] = current_value
        return results


    fig, ax = plt.subplots()

    # Floating-point configurations: (exp_bits, man_bits, bias, quantization)
    # Note: when the bias is set to None, it defaults to an IEEE-like bias
    fp_formats = [
        (8, 23, None, "TIES_EVEN"),
        (5, 10, None, "TIES_EVEN"),
        (5, 10, None, "TO_ZERO"),
        (5, 10, None, "STOCH_WEIGHTED"),
        (5, 10, None, "STOCH_EQUAL"),
        (4, 6, 11, "STOCH_WEIGHTED"),
    ]

    for exp_bits, man_bits, bias, mode in fp_formats:
        with apy.APyFloatQuantizationContext(
            eval(f"apy.QuantizationMode.{mode}"), quantization_seed=12
        ):
            results = solve_ode(exp_bits, man_bits, bias)
        ax.plot(results[:, 0], results[:, 1], label=f"E{exp_bits}M{man_bits} ({mode})")

    ax.set_aspect("equal")
    fig.legend(ncols=2, loc="upper center")
    plt.show()


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

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


.. _sphx_glr_download_examples_ode_solver.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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