Die Interpolation ist ein Teilgebiet der numerischen Mathematik und es befasst sich mit der Problematik, eine analytische Funktion mittels einer gegebenen Menge von Werten zu bestimmen. In der fünften Vorlesung hatten wir im Unterpunkt Anwendungsbeispiel: Interpolation und Polynomapproximation den wohl elementarsten Algorithmus einer Interpolation vorgestellt, bei dem ein Polynom $P(x)$ mittels einer gegebenen Liste von Werten konstruiert wird. Aufgrund des Weierstraßschen Approximationssatzes gibt es für jede, auf einem Teilintervall $[a,b]$ definierte, stetige Funktion, $f(x)$ ein Polynom $P(x)$ mit der Eigenschaft, dass für jedes $\epsilon > 0$ die Differenz der Funktion mit dem Polynom verschwindend klein wird ($\left| f(x) - P(x) \right| < \epsilon \,\, \forall \, x \in [a,b]$). Die Taylorschen Polynome sind hierbei wohl das bekannteste Beispiel einer Entwicklung einer Funktion in ein Polynom. Leider haben die Taylorschen Polynome den Nachteil, dass sie von ihrer Konstruktion her nur eine Stützstelle $x_0 \in [a,b]$ verwenden und somit die Funktion $f(x)$ zwar gut in dem Bereich um die Stützstelle beschreiben, aber nicht im gesamten Teilintervall $[a,b]$. Bei einer Interpolation mittels der Methode der Lagrange Polynome verwendet man hingegen mehrere Stützstellenpunkte für die Konstruktion des approximierenden Polynoms. Dieses Verfahren ist Gegenstand dieses Jupyter Notebooks.
Bei der Entwicklung einer Funktion in ein Lagrange Polynom steht man vor der Aufgabe eine Funktion $f(x)$ mittels ($n+1$) vorgegebener Punkte ($(x_0,f(x_0)), (x_1,f(x_1)), ..., (x_n,f(x_n))$) in ein Lagrange Polynom vom Grade $n$ zu entwickeln. Das n-te Lagrange Polynom $P_n(x)$ einer Funktion $f(x)$ ist wie folgt definiert: \begin{equation} P_n(x) = \sum_{k=0}^n f(x_k) \cdot L_{n,k}(x) \quad, \, \hbox{wobei} \quad L_{n,k}(x) = \prod_{i=0 , \, i \neq k}^n \frac{(x-x_i)}{(x_k-x_i)} \quad. \end{equation} Approximiert man z.B. die Funktion $f(x)=\frac{1}{x}$ im Teilintervall $[a,b]=[0.5,6]$ mittels dreier vorgegebener Punkte ($n=2$, mit $n+1$ Punkten $(2,f(2)), (2.5,f(2.5))$ und $(4,f(4)) $) in ein Lagrange Polynom vom Grade 2, so erhält man: \begin{eqnarray} P_2(x) = \sum_{k=0}^2 f(x_k) \cdot L_{2,k}(x) &=& f(x_0) \cdot L_{2,0}(x) + f(x_1) \cdot L_{2,1}(x) + f(x_2) \cdot L_{2,2}(x) \\ &=& \frac{1}{2} \cdot L_{2,0}(x) + \frac{1}{2.5} \cdot L_{2,1}(x) + \frac{1}{4} \cdot L_{2,2}(x) \\ \hbox{mit:} && L_{2,0}(x) = \prod_{i=0 , \, i \neq 0}^2 \frac{(x-x_i)}{(x_0-x_i)} = \frac{(x-x_1)}{(x_0-x_1)} \cdot \frac{(x-x_2)}{(x_0-x_2)} = \frac{(x-2.5)}{(2-2.5)} \cdot \frac{(x-4)}{(2-4)} = \frac{(x-2.5)}{(-0.5)} \cdot \frac{(x-4)}{(-2)} \\ && L_{2,1}(x) = \prod_{i=0 , \, i \neq 1}^2 \frac{(x-x_i)}{(x_1-x_i)} = \frac{(x-x_0)}{(x_1-x_0)} \cdot \frac{(x-x_2)}{(x_1-x_2)} = \, ... \\ && L_{2,2}(x) = \prod_{i=0 , \, i \neq 2}^2 \frac{(x-x_i)}{(x_2-x_i)} = \frac{(x-x_0)}{(x_2-x_0)} \cdot \frac{(x-x_1)}{(x_2-x_1)} = \, ... \\ \Longrightarrow P_2(x) = 0.05 \, x^2 - 0.425 \, x + 1.15 && \end{eqnarray}
In der Vorlesung hatten wir die Lagrange Polynom Methode in einem C++ Programm implementiert (siehe Lagrange_Polynom_1.cpp). Die Ergebnisse dieses Programms möchten wir zunächst hier visualisieren.
Die Terminalausgabe der Ergebnisse des C++ Programm wurde mittels des Befehls "./a.out > Lagrange_Polynom_1.dat" in eine Textdatei umgeleitet. Diese Textdatei lesen wir nun ein und speichern die zuerst ausgegebenen Stützstellenpunkte in einer Liste (Variablenname 'points'). Die eigentlichen Ergebnisdaten, welche erst nach der 7. Zeile beginnen, speichern wir unter dem Variablennamen 'data'.
Damit die folgenden Befehle ohne Fehlermeldung funktionieren, muss sich die Datei "Lagrange_Polynom_1.dat" in dem Jupyter Notebook Verzeichnis befinden.
import matplotlib.pyplot as plt
import numpy as np
points = np.genfromtxt("./Lagrange_Polynom_1.dat", max_rows=1)
data = np.genfromtxt("./Lagrange_Polynom_1.dat", skip_header=7)
Die untere Abbildung visualisiert die eingelesenen Daten des C++ Programms, wobei die blauen Punkte die drei vorgegebenen Stützstellen des Polynoms darstellen, bei welchen die Übereinstimmung $f(x_k)=P_2(x_k)$ exakt gilt.
plt.title(r'Entwicklung einer Funktion in ein Lagrange Polynom')
plt.ylabel(r'$f(x)\, , \,\, P_2(x)$')
plt.xlabel('x')
plt.plot(data[:,1],data[:,2], color="black", label=r'$f(x)$')
plt.plot(data[:,1],data[:,3], color="red", label=r'$P_2(x)$')
plt.scatter(points,1/points, marker='o', color="blue", s=25);
plt.legend(frameon=True, loc="upper right",fontsize=12);
In diesem Jupyter Notebook möchten wir außerdem den analytischen Ausdruck des Polynoms $P_2(x)$ mittels computerunterstützter symbolischer Berechnungen erstellen.
Zunächst wird das Python Modul 'sympy' eingebunden, das ein Computer-Algebra-System für Python bereitstellt und eine Vielzahl an symbolischen Berechnungen im Bereich der Mathematik und Physik relativ einfach möglich macht. Falls Sie das "sympy" Modul das erste Mal verwenden, müssen Sie es zunächst in Ihrer Python 3 Umgebung installieren (z.B. in einem Linux Terminal mit "pip3 install sympy").
from sympy import *
init_printing()
Wir möchten nun die Funktion $f(x)=\frac{1}{x}$ mittels ($n+1$) vorgegebener Punkte ($(x_0,f(x_0)), (x_1,f(x_1)), ..., (x_n,f(x_n)), $) in ein Lagrange Polynom vom Grade N entwickeln. Hierzu definieren wir zunächst das sympy-Symbol 'x' und die Funktion 'f':
x = Symbol('x')
f = 1/x
Das n-te Lagrange Polynom $P_n(x)$ einer Funktion ist wie folgt definiert:
$$ \begin{equation} P_n(x) = \sum_{k=0}^n f(x_k) \cdot L_{n,k}(x) \quad, \end{equation} $$mit
$$ \begin{equation} L_{n,k}(x) = \prod_{i=0 , \, i \neq k}^n \frac{(x-x_i)}{(x_k-x_i)} \quad. \end{equation} $$Wir möchten zunächst wieder die Funktion $f(x)=\frac{1}{x}$ mittels drei vorgegebener Punkte ($(2,f(2)), (2.5,f(2.5))$ und $(4,f(4)) $) in ein Lagrange Polynom vom Grade 2 approximieren. Wir definieren uns dafür die drei Stützstellenpunkte:
n=2
points=[2, 2.5, 4]
f_points=[f.subs(x,points[0]), f.subs(x,points[1]), f.subs(x,points[2])]
Den Algorithmus der Methode der Lagrange Polynome, welcher im C++ Programm implementiert wurde (siehe Lagrange_Polynom_1.cpp), kann man wie folgt in Python übertragen.
P_n = 0
for k in range(n+1):
Lk=1
for i in range(n+1):
if i != k:
Lk = Lk * (x - points[i])/(points[k] - points[i])
P_n = P_n + f_points[k] * Lk
Im Gegensatz zum C++ Programm, wird das approximierte Polynom $P_2(x)$ nun jedoch nicht als eine Zahlenliste ausgegeben, sondern mittels des Moduls 'sympy' als ein analytischer Ausdruck dargestellt. Es ergibt sich
P_n
und nach Vereinfachung des oberen Ausdrucks, die erwartete analytische Form für $P_2(x)$
P_n.simplify()
Man hätte den analytischen Ausdruck auch mittels der Python Funktionen 'sum(...)' und 'prod(...)' berechnen können:
List_ranges=[]
for i in range(n+1):
List_range=list(range(n+1))
List_range.pop(i)
List_ranges.append(List_range)
P_n=sum(f_points[k] * prod((x - points[i])/(points[k] - points[i]) for i in List_ranges[k]) for k in range(n+1))
P_n.simplify()
Wir stellen uns im Folgenden den analytischen Ausdruck von $P_2(x)$ zusammen mit der ursprünglichen Funktion $f(x)$ dar.
x_min=0.5
x_max=6
pf = plot(f, P_n.simplify(), (x,x_min,x_max), line_color='black', legend = True, show=False)
pf[1].line_color='red'
pf.show()
Man sieht, dass der gewonnene analytische Ausdruck identisch mit den berechneten Ergebnisse des C++ Programms ist.
Wir berechnen nun das 6-te Lagrange Polynom $P_{6}(x)$ der Funktion $f(x)=\frac{1}{x}$, wobei die sieben Stützstellen an den folgenden x-Werten genommen wurden: $\vec{x}=\left( 1, 1.5, 2, 2.5, 3, 5, 7 \right)$.
n=6
points=[1, 1.5, 2, 2.5, 3, 5, 7]
f_points=[]
for i in range(n+1):
f_points.append(f.subs(x,points[i]))
List_ranges=[]
for i in range(n+1):
List_range=list(range(n+1))
List_range.pop(i)
List_ranges.append(List_range)
P_n=sum(f_points[k] * prod((x - points[i])/(points[k] - points[i]) for i in List_ranges[k]) for k in range(n+1))
Nach Vereinfachung erhält man das folgende Polynom vom Grade $n=6$.
P_n.simplify()
x_min=0.5
x_max=7.5
pf = plot(f, P_n.simplify(), (x,x_min,x_max), line_color='black', show=False)
pf[1].line_color='red'
pf.show()
Wir vergleichen die zuvor dargestellten analytischen Berechnungen mit den simulierten Daten des C++ Programms ( Lagrange_Polynom_2.cpp ). Damit die folgenden Befehle ohne Fehlermeldung funktionieren, muss sich die Datei "Lagrange_Polynom_2.dat" in dem Jupyter Notebook Verzeichnis befinden.
points = np.genfromtxt("./Lagrange_Polynom_2.dat", max_rows=1)
data = np.genfromtxt("./Lagrange_Polynom_2.dat", skip_header=7)
plt.title(r'Entwicklung einer Funktion in ein Lagrange Polynom')
plt.ylabel(r'$f(x)\, , \,\, P_6(x)$')
plt.xlabel('x')
plt.plot(data[:,1],data[:,2], color="black", label=r'$f(x)$')
plt.plot(data[:,1],data[:,3], color="red", label=r'$P_6(x)$')
plt.scatter(points,1/points, marker='o', color="blue", s=25);
plt.legend(frameon=True, loc="upper right",fontsize=12);
Man sieht wieder, dass die berechneten Ergebnisse identisch sind.