[NumCS] Add code (adapt. quad)

semester3/numcs/parts/02_quadrature/03_adaptive.tex

@@ -3,6 +3,7 @@
 % └                                                ┘
 
 \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*}
@@ -17,3 +18,31 @@ Auf Seiten 150 - 151 im Skript findet sich Code, um eine adaptive Quadratur durc
 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}

semester3/numcs/parts/02_quadrature/04_in-rd.tex

@@ -23,8 +23,6 @@ Was dasselbe ist, wie oben, aber mit $d$ Summen und $d$-mal ein $w_{j_k}$ und ei
 
 % https://www.slingacademy.com/article/scipy-integrate-simpson-function-4-examples/ explains scipy's n-d integration well
 
-\newpage
-
 \innumpy Lassen sich so die bekannten Verfahren wie die Trapez-Regel oder Simpson-Regel leicht auf höhere Dimensionen anwenden
 
 \begin{code}{python}
@@ -64,6 +62,8 @@ def simpson_2d_mesh(f, a: float, b: float, Nx: int, c: float, d: float, Ny: int)
     return Q
 \end{code}
 
+\newpage
+
 \begin{recall}[]{Tensor-Produkt}
     \TODO Write this section
 \end{recall}