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


Simpson's Newton method for systems of equations
================================================

Thomas Simpson's 1740 *Essays* stated Newton's method in terms of
fluxions and applied it to two equations in two unknowns. In matrix form,
linearize :math:`F: \mathbb R^n \to \mathbb R^n` at the current guess and
solve for the step:

.. math::

   J(\mathbf x_k)\,\mathbf s_k = -F(\mathbf x_k), \qquad
   \mathbf x_{k+1} = \mathbf x_k + \mathbf s_k .

Here :math:`F(x, y) = (x^2 + y^2 - 4,\; e^x + y - 1)`: the circle of
radius 2 meets the curve :math:`y = 1 - e^x` twice. Each starting point
is drawn to one intersection, and near it the residual is roughly squared
at every step.

.. GENERATED FROM PYTHON SOURCE LINES 22-38

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from mathematicskit.numerical_analysis import NewtonSystem


    def F(v):
        x, y = v
        return np.array([x**2 + y**2 - 4.0, np.exp(x) + y - 1.0])


    def J(v):
        x, y = v
        return np.array([[2.0 * x, 2.0 * y], [np.exp(x), 1.0]])









.. GENERATED FROM PYTHON SOURCE LINES 39-41

Iterate paths in the plane
--------------------------

.. GENERATED FROM PYTHON SOURCE LINES 41-60

.. code-block:: Python


    starts = [(-2.5, 2.5), (-0.5, 2.0), (2.5, 0.5), (0.5, -2.5), (-2.0, -1.0)]
    fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(12, 5))
    xs = np.linspace(-3, 3, 400)
    X, Y = np.meshgrid(xs, xs)
    ax0.contour(X, Y, X**2 + Y**2 - 4.0, levels=[0.0], colors="tab:blue")
    ax0.contour(X, Y, np.exp(X) + Y - 1.0, levels=[0.0], colors="tab:orange")
    for k, start in enumerate(starts):
        result = NewtonSystem(F, start, jacobian=J, tol=1e-14).solve()
        ax0.plot(result.history[:, 0], result.history[:, 1], "o-", ms=4, color=f"C{k + 2}")
        ax1.semilogy(result.residual_norms + 1e-17, "o-", color=f"C{k + 2}", label=f"start {start}")
        print(f"start {start}: root ({result.root[0]:+.12f}, {result.root[1]:+.12f}) in {result.iterations} iterations")
    ax0.set_xlim(-3, 3)
    ax0.set_ylim(-3, 3)
    ax0.set_aspect("equal")
    ax0.set_xlabel("x")
    ax0.set_ylabel("y")
    ax0.set_title(r"$x^2 + y^2 = 4$ (blue) meets $y = 1 - e^x$ (orange)")




.. image-sg:: /api/gallery/numerical_analysis/root_finding/images/sphx_glr_plot_05_newton_systems_001.png
   :alt: $x^2 + y^2 = 4$ (blue) meets $y = 1 - e^x$ (orange)
   :srcset: /api/gallery/numerical_analysis/root_finding/images/sphx_glr_plot_05_newton_systems_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none

    start (-2.5, 2.5): root (-1.816264068825, +0.837367799891) in 7 iterations
    start (-0.5, 2.0): root (-1.816264068825, +0.837367799891) in 7 iterations
    start (2.5, 0.5): root (-1.816264068825, +0.837367799891) in 9 iterations
    start (0.5, -2.5): root (+1.004168738475, -1.729637287026) in 6 iterations
    start (-2.0, -1.0): root (-1.816264068825, +0.837367799891) in 7 iterations

    Text(0.5, 1.0, '$x^2 + y^2 = 4$ (blue) meets $y = 1 - e^x$ (orange)')



.. GENERATED FROM PYTHON SOURCE LINES 61-63

Quadratic convergence of the residual
-------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 63-71

.. code-block:: Python


    ax1.set_xlabel("iteration k")
    ax1.set_ylabel(r"$\|F(\mathbf{x}_k)\|_2$")
    ax1.set_title("Correct digits roughly double per step")
    ax1.legend(fontsize=8)
    fig.tight_layout()

    plt.show()








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

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


.. _sphx_glr_download_api_gallery_numerical_analysis_root_finding_plot_05_newton_systems.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/numerical_analysis/root_finding/plot_05_newton_systems.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

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

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

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

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

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

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


.. only:: html

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

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