mirror of
https://github.com/janishutz/eth-summaries.git
synced 2026-03-14 23:10:03 +01:00
48 lines
2.3 KiB
TeX
48 lines
2.3 KiB
TeX
% ┌ ┐
|
|
% │ 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} |