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


The Thomas algorithm
====================

Every implicit step in one space dimension is a tridiagonal system.
Llewellyn Thomas (1949) wrote down the specialised Gaussian elimination
that solves it in :math:`O(m)` operations: one sweep removes the
sub-diagonal, a second substitutes back. No pivoting is needed because
the implicit Black-Scholes matrix is diagonally dominant.
:func:`~financialkit.pde_methods.thomas_solve` is that algorithm by hand;
it agrees with :func:`scipy.linalg.solve_banded` to rounding, and its cost
grows linearly while a dense solve grows cubically.

.. GENERATED FROM PYTHON SOURCE LINES 16-34

.. code-block:: Python

    import time

    import matplotlib.pyplot as plt
    import numpy as np

    from financialkit.pde_methods import CrankNicolson, black_scholes_operator, solve_tridiagonal, thomas_solve, tridiagonal_matrix

    # One implicit step's matrix, I - dt A, on a 200-interval grid
    S = np.linspace(0.0, 400.0, 201)
    lower, diag, upper = black_scholes_operator(S, 0.2, 0.05)
    dt = 0.01
    lo, di, up = -dt * lower, 1.0 - dt * diag, -dt * upper
    rhs = np.maximum(S[1:-1] - 100.0, 0.0)

    x = thomas_solve(lo, di, up, rhs)
    print("max |Thomas - solve_banded| =", np.abs(x - solve_tridiagonal(lo, di, up, rhs)).max())
    print("diagonally dominant:", bool(np.all(np.abs(di) >= np.abs(lo) + np.abs(up))))





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

 .. code-block:: none

    max |Thomas - solve_banded| = 1.7053025658242404e-13
    diagonally dominant: True




.. GENERATED FROM PYTHON SOURCE LINES 35-37

Linear against cubic cost
-------------------------

.. GENERATED FROM PYTHON SOURCE LINES 37-58

.. code-block:: Python


    sizes = np.array([100, 200, 400, 800, 1600])
    t_thomas, t_dense = [], []
    for m in sizes:
        lo_m, up_m = np.full(m, -1.0), np.full(m, -1.0)
        di_m, b = np.full(m, 2.5), np.ones(m)
        dense = tridiagonal_matrix(lo_m, di_m, up_m).toarray()
        t0 = time.perf_counter()
        thomas_solve(lo_m, di_m, up_m, b)
        t_thomas.append(time.perf_counter() - t0)
        t0 = time.perf_counter()
        np.linalg.solve(dense, b)
        t_dense.append(time.perf_counter() - t0)

    fig, axes = plt.subplots(1, 2, figsize=(11, 4))
    axes[0].loglog(sizes, t_thomas, "o-", label="Thomas, O(m)")
    axes[0].loglog(sizes, t_dense, "s-", label="dense LU, O(m^3)")
    axes[0].set_xlabel("unknowns m")
    axes[0].set_ylabel("seconds")
    axes[0].legend()




.. image-sg:: /api/gallery/pde_methods/schemes/images/sphx_glr_plot_03_thomas_algorithm_001.png
   :alt: plot 03 thomas algorithm
   :srcset: /api/gallery/pde_methods/schemes/images/sphx_glr_plot_03_thomas_algorithm_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none


    <matplotlib.legend.Legend object at 0x116c09e80>



.. GENERATED FROM PYTHON SOURCE LINES 59-61

The same prices from either solver
----------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 61-69

.. code-block:: Python


    a = CrankNicolson(200, 50, solver="thomas").price(100.0, 100.0, 1.0, 0.05, 0.2)
    b = CrankNicolson(200, 50, solver="banded").price(100.0, 100.0, 1.0, 0.05, 0.2)
    axes[1].semilogy(a.S, np.abs(a.V - b.V) + 1e-18)
    axes[1].set_xlabel("spot S")
    axes[1].set_title("|Thomas - solve_banded| over the grid")
    fig.tight_layout()
    plt.show()








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

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


.. _sphx_glr_download_api_gallery_pde_methods_schemes_plot_03_thomas_algorithm.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/pde_methods/schemes/plot_03_thomas_algorithm.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

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

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

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

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

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

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


.. only:: html

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

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