% ┌                                                ┐
% │     AUTHOR: Janis Hutz<info@janishutz.com>     │
% └                                                ┘

\newsection

\subsection{Adaptive Quadratur}
Der lokale Fehler einer zusammengesetzten Quadraturformel auf dem Gitter $\mathcal{M} := \{ a = x_0 < x_1 < \dots < x_m = b \}$ ist (für $f \in C^2([a, b])$):
\begin{align*}
    \left| \int_{x_k}^{x_{k + 1}} f(t) \dx t - \frac{f(x_k) + f(x_{k + 1})}{2}(x_{k + 1} - x_k) \right| \leq (x_{k + 1} - x_k)^3 ||f''||_{L^\infty([x_k, x_{k + 1}])}
\end{align*}
Also ist es nur sinnvoll, das Gitter zu verfeinern wo $|f''|$ gross ist.

Auf Seiten 150 - 151 im Skript findet sich Code, um eine adaptive Quadratur durchzuführen.

\setLabelNumber{all}{3}
\fancyremark{Adaptive Quadratur in Python} Mit \texttt{scipy.integrate.quad} können wir einfach eine adaptive Quadratur durchführen und benutzt \texttt{QUADPACK}.
Mit \texttt{scipy.integrate.quadrature} können wir die Gauss-Quadratur verwenden.

Für $x \in \R^d$, also eine mehrdimensionale Funktion der Dimension $d$ können wir \texttt{scipy.integrate.nquad} verwenden. Mehr dazu im nächsten Kapitel

\innumpy Lässt sich eine simple adaptive Quadratur folgendermassen implementieren. Die Idee ist hierbei einfach das Gitter an Problemstellen zu verdichten, was iterativ die Konvergenz der Quadratur verbessern sollte.

\begin{code}{python}
def adaptive_quad(f, M: np.ndarray, rtol: float = 1e-6, atol: float = 1e-10):
    """ Calculate adaptive quad. of f on given grid M"""
    h = np.diff(M)
    midp = 0.5 * ( M[:-1] + M[1:] )
    fx   = f(M); fmid = f(midp)
    
    # Local quadratures to estimate local errors further down
    trap_local = h/4 * ( fx[:-1] + 2*fmid + fx[1:] )
    simp_local = h/6 * ( fx[:-1] + 4*fmid + fx[1:] )
    Q = np.sum(simp_local)
    
    # Estimate local & total error
    err_loc = np.abs( simp_local - trap_local )
    err = np.sum(err_loc)
    err_avg = err / len(err_loc)
    
    # Refine grid in high-error regions, if needed
    if (err > rtol*np.abs(Q) and err > atol):
        refcells = np.nonzero( err_loc > 0.9 * err_avg )[0]     # Find problematic cells
        M_new = np.sort( np.append(M, midp[refcells]) )         # Update grid
        Q = adaptive_quad(f, M_new, rtol, atol, results)        # Try again
    
    return Q
\end{code}