diff --git a/semester3/numcs/numcs-summary.pdf b/semester3/numcs/numcs-summary.pdf index 97a5e58..2232fda 100644 Binary files a/semester3/numcs/numcs-summary.pdf and b/semester3/numcs/numcs-summary.pdf differ diff --git a/semester3/numcs/parts/02_quadrature/03_adaptive.tex b/semester3/numcs/parts/02_quadrature/03_adaptive.tex index b183cef..4ddd478 100644 --- a/semester3/numcs/parts/02_quadrature/03_adaptive.tex +++ b/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} \ No newline at end of file diff --git a/semester3/numcs/parts/02_quadrature/04_in-rd.tex b/semester3/numcs/parts/02_quadrature/04_in-rd.tex index 225fd2f..769c71e 100644 --- a/semester3/numcs/parts/02_quadrature/04_in-rd.tex +++ b/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}