.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/linalg/lu/plot_01_nine_chapters_elimination.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_linalg_lu_plot_01_nine_chapters_elimination.py>`
        to download the full example code.

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

.. _sphx_glr_api_gallery_linalg_lu_plot_01_nine_chapters_elimination.py:


Elimination in the Nine Chapters: the three grades of grain
===========================================================

The first problem of Chapter 8 (*Fangcheng*, "rectangular arrays") of
the *Nine Chapters on the Mathematical Art*: 3 bundles of top-grade,
2 of medium-grade and 1 of low-grade grain yield 39 *dou*; 2, 3 and 1
bundles yield 34; 1, 2 and 3 bundles yield 26. How much does one bundle
of each grade yield?

The text lays the coefficients out on a counting board and eliminates
unknowns without ever dividing: to clear an entry it multiplies a whole
row by the pivot and subtracts a multiple of the pivot row. This script
replays that division-free elimination step by step and checks the
answer against a modern pivoted LU solve.

.. GENERATED FROM PYTHON SOURCE LINES 19-25

.. code-block:: Python

    import matplotlib.pyplot as plt
    import numpy as np

    from mathematicskit.linalg import lu_solve_system
    from mathematicskit.linalg.visualizers.plots import plot_matrix_heatmap








.. GENERATED FROM PYTHON SOURCE LINES 26-29

The counting board
------------------
One row per bundle mix: [top, medium, low | yield].

.. GENERATED FROM PYTHON SOURCE LINES 29-39

.. code-block:: Python


    board = np.array(
        [
            [3.0, 2.0, 1.0, 39.0],
            [2.0, 3.0, 1.0, 34.0],
            [1.0, 2.0, 3.0, 26.0],
        ]
    )
    print("initial board:\n", board.astype(int))





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

 .. code-block:: none

    initial board:
     [[ 3  2  1 39]
     [ 2  3  1 34]
     [ 1  2  3 26]]




.. GENERATED FROM PYTHON SOURCE LINES 40-44

Division-free elimination, column by column
-------------------------------------------
To clear ``board[j, i]`` the row ``j`` is replaced by
``pivot * row_j - board[j, i] * row_i`` -- integers stay integers.

.. GENERATED FROM PYTHON SOURCE LINES 44-53

.. code-block:: Python


    stages = [board.copy()]
    for i in range(2):
        pivot = board[i, i]
        for j in range(i + 1, 3):
            board[j] = pivot * board[j] - board[j, i] * board[i]
        stages.append(board.copy())
        print(f"after clearing column {i}:\n", board.astype(int))





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

 .. code-block:: none

    after clearing column 0:
     [[ 3  2  1 39]
     [ 0  5  1 24]
     [ 0  4  8 39]]
    after clearing column 1:
     [[ 3  2  1 39]
     [ 0  5  1 24]
     [ 0  0 36 99]]




.. GENERATED FROM PYTHON SOURCE LINES 54-56

Back substitution
-----------------

.. GENERATED FROM PYTHON SOURCE LINES 56-67

.. code-block:: Python


    x = np.zeros(3)
    for i in (2, 1, 0):
        x[i] = (board[i, 3] - board[i, i + 1 : 3] @ x[i + 1 :]) / board[i, i]
    print("yield per bundle (top, medium, low):", x, "dou")
    print("Nine Chapters answer: 9 1/4, 4 1/4, 2 3/4 dou")

    A = stages[0][:, :3]
    b = stages[0][:, 3]
    print("agrees with pivoted LU solve:", np.allclose(x, lu_solve_system(A, b)))





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

 .. code-block:: none

    yield per bundle (top, medium, low): [9.25 4.25 2.75] dou
    Nine Chapters answer: 9 1/4, 4 1/4, 2 3/4 dou
    agrees with pivoted LU solve: True




.. GENERATED FROM PYTHON SOURCE LINES 68-70

The board becomes triangular
----------------------------

.. GENERATED FROM PYTHON SOURCE LINES 70-82

.. code-block:: Python


    fig, axes = plt.subplots(1, 3, figsize=(10, 3.2))
    titles = ["initial board", "column 0 cleared", "column 1 cleared (triangular)"]
    for ax, stage, title in zip(axes, stages, titles):
        plot_matrix_heatmap(stage, ax=ax, title=title)
        for (r, c), val in np.ndenumerate(stage):
            ax.text(c, r, f"{val:.0f}", ha="center", va="center", fontsize=9)
        ax.set_xticks(range(4), ["top", "med", "low", "yield"])
        ax.set_yticks(range(3))
    fig.tight_layout()

    plt.show()



.. image-sg:: /api/gallery/linalg/lu/images/sphx_glr_plot_01_nine_chapters_elimination_001.png
   :alt: initial board, column 0 cleared, column 1 cleared (triangular)
   :srcset: /api/gallery/linalg/lu/images/sphx_glr_plot_01_nine_chapters_elimination_001.png
   :class: sphx-glr-single-img






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

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


.. _sphx_glr_download_api_gallery_linalg_lu_plot_01_nine_chapters_elimination.py:

.. only:: html

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

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

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

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

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

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

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


.. only:: html

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

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