[NumCS] Update/Add code examples

semester3/numcs/parts/01_interpolation/00_polynomial/01_monome.tex

@@ -30,10 +30,28 @@ $$
 	\end{bmatrix}
 $$
 
+\begin{code}{python}
+def coeffs_monomial(x: np.ndarray, y: np.ndarray):
+    """ Solve Vandermonde matrix for monomial coeffs (very unstable) """
+    A = np.vander(x)
+    coeffs = np.linalg.solve(A, y)
+    return coeffs
+\end{code}
+
 Um $\alpha_i$ zu finden ist die Vandermonde Matrix unbrauchbar, da die Matrix schlecht konditioniert ist.
 
+\newpage
+
 Zur Auswertung von $p(x)$ kann man direkt die Matrix-darstellung nutzen, oder effizienter:
 
 \fancydef{Horner Schema} $p(x) = (x \ldots x ( x (\alpha_n x + \alpha_{n-1}) + \ldots + \alpha_1) + \alpha_0)$
 
+\begin{code}{python}
+def eval_horner(coeffs: np.ndarray, vals: np.ndarray):
+    """ Evaluate polynomial using Horner scheme """
+    h = coeffs[0]
+    for i in range(1, len(coeffs)): h = vals * h + coeffs[i]
+    return h
+\end{code}
+
 \innumpy liefert \verb|polyfit| die direkte Auswertung, \verb|polyval| wertet Polynome via Horner-Schema aus. (Gemäss Script, in der Praxis sind diese Funktionen \verb|deprecated|)

semester3/numcs/parts/01_interpolation/00_polynomial/02_newton-basis.tex

@@ -1,4 +1,3 @@
-\newpage
 \subsection{Newton Basis}
 % Session: Herleitung unwichtig, konzentrieren auf Funktion/Eigenschaften von Newton/Lagrange.
 
@@ -91,8 +90,6 @@ Falls $x_j = x_0 + \underbrace{j \cdot h}_{:= \Delta^j}$ gilt vereinfacht sich e
 wenn bspw. die Stützstellen nicht gut gewählt sind oder das Polynom einen zu hohen Grad hat.
 Sie ist definiert durch $\displaystyle f(x) = \frac{1}{1 + x^2}$
 
-
-\newpage
 \begin{multicols}{2}
 
     Matrixmultiplikation in $\mathcal{O}(n^3)$, Speicher $\mathcal{O}(n^2)$
@@ -171,7 +168,3 @@ Auswertung eines Newton-Polynoms funktioniert in $\mathcal{O}(n)$ durch ein modi
 $$
     \forall x \in [a,b]\ \exists \xi \in (a,b):\quad\quad \underbrace{f(x)-p(x)}_{\text{Fehler}} = \prod_{i=0}^{n}(x-x_i)\cdot\frac{f^{(n+1)}(\xi)}{(n+1)!}
 $$
-
-Man bemerke: Die Wahl der Stützpunkte hat direkten Einfluss auf den Fehler.
-
-

semester3/numcs/parts/01_interpolation/00_polynomial/03_lagrange-and-barzycentric-formula.tex

@@ -51,37 +51,15 @@ Man berechnet die baryzentrischen Gewichte $\lambda_k$ folgendermassen:
 \end{align*}
 oder das ganze mithilfe von Numpy:
 \begin{code}{python}
-    def barycentric_weights(x: np.ndarray) -> np.ndarray:
+def weights_barycentric(x: np.ndarray):
-        n = len(x)
+    """ All x should be pairwise distinct, else zero-divison error. """
-        # Initialize to zeros
+    n = x.size
-        barweight = np.ones(n)
+    w = np.zeros(n)
-        for k in range(n):
+    for i in range(n):
-            # Vectorized differences between $x_k$ and all $x$s
+        w[i] = 1/( np.prod(x[i] - x[0:i]) * np.prod(x[i] - x[i+1:n]) )
-            differences = x[k] - x
+    return w
-            # Remove the $k$-th element (and handle edge cases for $k = 0$ and $k = n - 1$)
-            if k < n - 1 and k > 0:
-                diff_processed = np.concatenate((differences[:k], differences[(k + 1) :]))
-                barweight[k] = 1 / np.prod(diff_processed)
-            elif k == 0:
-                barweight[k] = 1 / np.prod(differences[1:])
-            else:
-                barweight[k] = 1 / np.prod(differences[:k])
-        return barweight
 \end{code}
 
-Gleiche funktion, etwas kürzer:
-
-\begin{code}{python}
-    def barycentric_weights(x: np.ndarray) -> np.ndarray:
-        n = len(x)
-        w = np.ones(n) # = barweight
-        # Compute the (non-inverted) product, avoiding case (x[i] - x[i]) = 0
-        for i in range(0, n, 1):
-            if (i-1 > 0): w[0:(i-1)] *= (x[0:(i-1)] - x[i])
-            if (i+1 < n): w[i+1:n]   *= (x[i+1:n]   - x[i])
-        # Invert all at once
-        return 1/w
-\end{code}
 
 Mit dem können wir dann ein Polynom mit der baryzentrischen Interpolationsformel interpolieren:
 \numberingOff
@@ -101,26 +79,20 @@ Mit anderen $\lambda_k$ eröffnet die baryzentrische Formel einen Weg zur Verall
 Eine weitere Anwendung der Formel ist als Ausganspunkt für die Spektralmethode für Differenzialgleichungen.
 
 \begin{code}{python}
-    def interp_barycentric(
+def eval_barycentric(w: np.ndarray, data: np.ndarray, y: np.ndarray, x: np.ndarray):
-        data_point_x: np.ndarray,
+    """  Sequentially calculate the barycentric formula """
-        data_point_y: np.ndarray,
+    n = x.size
-        barweight: np.ndarray,
+    tmp = np.ones(n)
-        x: np.ndarray
+    for i in range(n): tmp[i] = eval_barycentric_scalar(w, data, y, x[i])
-    ):
+    return tmp
-        p_x = np.zeros_like(x)
-        n = data_point_x.shape[0]
     
-        for i in range(x.shape[0]):
-            # Separate sums to divide in the end
-            upper_sum = 0
-            lower_sum = 0
-            for k in range(n):
-                frac = barweight[k] / (x[i] - data_point_x[k])
-                upper_sum += frac * data_point_y[k]
-                lower_sum += frac
-            p_x[i] = upper_sum / lower_sum
     
-        return p_x
+def eval_barycentric_scalar(w: np.ndarray, data: np.ndarray, y: np.ndarray, x):
+    """ Barycentric interpolation formula for a single value x """
+    n = data.size;
+    bottom = np.sum( w / (x - data) )
+    top = np.sum( w / (x - data) * y)
+    return top / bottom
 \end{code}
 
 

semester3/numcs/parts/01_interpolation/00_polynomial/04_chebyshev-interpolation.tex

@@ -139,20 +139,15 @@ Auf der nächsten Seite findet sich eine saubere, effiziente Implementation des
 
 \newpage
 \begin{code}{python}
-    def clenshaw(coeffs: np.ndarray, x: np.ndarray):
+def clenshaw(c: np.ndarray, x: np.ndarray):
-        n = len(coeffs) - 1
+    """ Clenshaw algorithm to evaluate polynomial using chebyshev coeffs """
-        # initialise temporary variables
+    n = c.size; m = x.size
-        d_prev_prev, d_prev, d_curr = (
+    d = np.zeros((m, 3))    # Save vectors [curr, prev1, prev2] as matrix
-            np.zeros_like(x),
-            np.zeros_like(x),
-            np.zeros_like(x),
-        )
     
-        for k in range(n, -1, -1):  # backward recursion
+    for i in range(n-1, -1, -1):
-            d_curr = coeffs[k] + 2 * x * d_prev - d_prev_prev
+        d[:, 2] = d[:, 1]; d[:, 1] = d[:, 0]
-            d_prev_prev, d_prev = d_prev, d_curr
+        d[:, 0] = c[i] + (2*x)*d[:, 1] - d[:, 2]
-
+    return d[:, 0] - x*d[:, 1]
-        return d_prev - x * d_prev_prev
 \end{code}
 
 

semester3/numcs/parts/01_interpolation/01_trigonometric/03_interpolation/01_zero-padding.tex

@@ -8,21 +8,49 @@ Ein trigonometrisches Polynom $p_{N - 1}(t)$ kann effizient an den äquidistante
 Dazu muss das Polynom $p_{N - 1} \in \mathcal{T}_N \subseteq \mathcal{T}_M$ in der trigonometrischen Basis $\mathcal{T}_M$ neugeschrieben werden,
 in dem man \bi{Zero-Padding} verwendet, also Nullen im Koeffizientenvektor an den Stellen höheren Frequenzen einfügt.
 
-\TODO Insert cleaned up code from Page 95 (part of exercises)
+\innumpy Dieses Verfahren lässt sich in Python leicht via \verb|slices| umsetzen:
-
-Die folgende Funktion wird im Script \texttt{evaliptrig} genannt.
 \rmvspace
 \begin{code}{python}
-    def evaluate_trigonometric_interpolation_polynomial(y: np.ndarray, N: int):
+def zero_pad(v: np.ndarray, N: int):
-        n = len(y)
+    """Apply zero-padding to size N to a vector v """
-        if (n % 2) == 0:
+    n = v.size
-            c = np.fft.ifft(y)  # Fourier coefficients
+    if (N < n): raise ValueError(f"ERROR: Zeropadding for N smaller than vector length: {N} < {n}")
-            a = np.zeros(N, dtype=complex)
     
-            # Zero padding
+    u = np.zeros(N, dtype=complex)
-            a[: n // 2] = c[: n // 2]
+    u[:n//2] = v[:n//2]
-            a[N - n // 2 :] = c[n // 2 :]
+    u[N-n//2:] = v[n//2:]
-            return np.fft.fft(a)
+    return u
-        else:
+
-            raise ValueError("odd length")
+
+def eval_trig_poly(y: np.ndarray, N: int):
+    """ Evaluate trig poly generated using y on N points """
+    n = y.size
+    if (n % 2 != 0): raise ValueError(f"ERROR: y must be of even length, len(y)={n}")
+    
+    coeffs = np.fft.fft(y) * 1/n
+    coeffs = zero_pad(coeffs, N)
+    return np.fft.ifft(coeffs) * N
 \end{code}
+
+\subsubsection{Ableitungen mit Zero-Padding}
+Mit dem Trick aus Bemerkung 3.2.13 lassen sich auch direkt die Ableitungen berechnen.
+
+\innumpy Geht dies direkt durch leichte Modifikation der obigen Funktionen:
+
+\begin{code}{python}
+def eval_trig_poly_d1(y: np.ndarray, N: int):
+    """ Evaluates first der. of trig poly generated using y on N points """
+    n = y.size
+    if (n % 2 != 0): raise ValueError(f"ERROR: y must be of even length, len(y)={n}")
+    
+    coeffs = np.fft.fft(y) * 1/n
+    
+    for i in range(0, n//2):
+        coeffs[i] *= (2.0j * np.pi * i)
+    for i in range(n//2, n):
+        coeffs[i] *= (2.0j * np.pi * (i - n))
+    
+    coeffs = zero_pad(coeffs, N)
+    return np.fft.ifft(coeffs) * N
+\end{code}
+

semester3/numcs/parts/01_interpolation/02_piece-wise/00_intro.tex

@@ -39,3 +39,28 @@ also gilt:
 \end{align*}
 Dies ist eine lokale Interpolation und $s_j$ ist $0$ ausser im definierten Intervall. 
 Die Idee des Variablenwechsel ist es, das Intervall, auf welchem die Funktion definiert ist von $[0, 1]$ nach $[x_{j - 1}, x_j]$ zu verschieben.
+
+\innumpy Stückweise lineare interpolation lässt sich leicht umsetzen via numpy \verb|piecewise| Funktionen:
+
+\begin{code}{python}
+def slope(x: np.ndarray, j: int, x_data: np.ndarray, y: np.ndarray):
+    """ Slope on j-th sub-interval, for x falling inside that sub-interval """
+    h = np.abs(x_data[j] - x_data[j-1])   # size of sub-interval
+    return y[j-1]*(x_data[j] - x)/h + y[j]*(x - x_data[j-1])/h
+
+
+def eval_linear_interp(x: np.ndarray, x_data: np.ndarray, y: np.ndarray):
+    """ Evaluate linear interpolation on x, using data points (x_data, y) """
+    n = x_data.size; m = x.size
+        
+    condlist = [
+        (x_data[j-1] <= x) & (x <= x_data[j])
+        for j in range(1, n)
+    ]
+    funclist = [
+        lambda x, ind=j: slope(x, ind, x_data, y)  
+        for j in range(1, n)
+    ]
+    
+    return np.piecewise(x, condlist, funclist)
+\end{code}