Compare a few numerical integration techniques and show off chebyshev-gauss quadrature

chebyshev_gauss.py

import marimo

__generated_with = "0.13.13-dev18"
app = marimo.App()


@app.cell
def _():
    import marimo as mo
    import numpy as np
    import matplotlib.pyplot as plt
    from scipy import integrate
    return integrate, mo, np, plt


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        """
        **This  notebook is powered by [WASM](https://webassembly.org/)**:
        it's running entirely in your browser!

        """
    ).callout()
    return


@app.cell(hide_code=True)
def _(mo):
    mo.md(
        r"""
    # chebyshev polynomial roots
    Note that they are concentrated on the edges and in the interval [-1,1].
    """
    )
    return


@app.cell(hide_code=True)
def _(mo):
    slider = mo.ui.slider(10, 40, show_value=True, label="Number of nodes")
    slider
    return (slider,)


@app.cell(hide_code=True)
def _(np, plt, slider):
    n = slider.value
    nodes = np.cos(np.pi * (0.5 + np.arange(0, n)) / n)
    # note it's between [-1,1]
    plt.eventplot(nodes)
    plt.title(f"Chebyshev nodes, n={n}")
    return (n,)


@app.cell(hide_code=True)
def _(mo):
    mo.md(r"""# numerical approximation of integral""")
    return


@app.cell(hide_code=True)
def _(ans, err_rate, mo, n):
    mo.md(
        f"""
    |     | Value | Error %
    | -------- | ------- | ------- |
    | *Ground Truth*  |  {ans['ground truth']:.5f}   | 0 |
    | Chebyshev-Trapezoid {n} nodes | {ans['approximation: chebyshev trapezoid']:.5f} | {err_rate['approximation: chebyshev trapezoid']:.5f}%| 
    | Chebyshev-Gauss {n} nodes  | {ans['approximation: chebyshev gauss']:.5f}|{err_rate['approximation: chebyshev gauss']:.5f}%|

    """
    )
    return


@app.cell
def _(integrate, n, np):
    def integrate_chebyshev_gauss(func, a, b):
        # [-1,1]
        chebyshev_nodes = np.cos(np.pi * (0.5 + np.arange(0, n)) / n)
        # [a,b]
        y = ((b - a) / 2) * chebyshev_nodes + ((b - a) / 2) + a
        # this is actually the jacobian!
        dy_coeff = (b - a) / 2  # times dx (x is chebyshev nodes)
        chebyshev_func_form_coeff = np.sqrt(1 - chebyshev_nodes**2)
        weights = np.pi / n

        return np.sum(weights * func(y) * dy_coeff * chebyshev_func_form_coeff)


    def integral_chebyshev_trapezoid(func, a, b):
        # x = np.linspace(a,b,num=n)
        chebyshev_nodes = np.cos(np.pi * (0.5 + np.arange(0, n)) / n)
        x = ((b - a) / 2) * chebyshev_nodes + ((b - a) / 2) + a
        x = sorted(x) # trapezoid integration needs ordered points
        return integrate.trapezoid(func(x), x=x)

    def integral_groundtruth(a, b):
        return -np.cos(b) - (-np.cos(a))


    a, b = 0, np.pi
    func = np.sin

    ans_gt = integral_groundtruth(a,b)
    ans_tr = integral_chebyshev_trapezoid(func, a, b)
    ans_cg = integrate_chebyshev_gauss(func, a, b)

    ans = {
        'ground truth':ans_gt,
        'approximation: chebyshev trapezoid':ans_tr,
        'approximation: chebyshev gauss':ans_cg,
    }

    err_rate = {
        'ground truth':0,
        'approximation: chebyshev trapezoid':(ans_gt-ans_tr)*100/ans_gt,
        'approximation: chebyshev gauss':(ans_gt-ans_cg)*100/ans_gt, 
    }

    ans

    return ans, err_rate


@app.cell
def _(ans, n, plt):
    labels = list(ans.keys())
    values = list(ans.values())

    plt.figure(figsize=(20, 6))
    plt.bar(labels, values, color=['skyblue', 'lightcoral', 'lightgreen'])
    plt.yscale("log")
    plt.ylabel("Value (log scale)")
    plt.title(f"Integral approximations with various methods - using {n} nodes")
    plt.show()
    return


@app.cell
def _():
    return


if __name__ == "__main__":
    app.run()