Commit 35bfebd4 authored by Christoph Ruegge's avatar Christoph Ruegge
Browse files

Merge remote-tracking branch 'refs/remotes/origin/master'

parents 72a47df0 e9a13a54
......@@ -40,7 +40,6 @@ Folgende Themen werden vermittelt:
* [18_Nichtlineare Gleichungen](http://nbviewer.jupyter.org/urls/gitlab.gwdg.de/jschulz1/math_prog/raw/master/lecture/18_Nichtlineare Gleichungen.ipynb)
* [19_Objektorientierung](http://nbviewer.jupyter.org/urls/gitlab.gwdg.de/jschulz1/math_prog/raw/master/lecture/19_Objektorientierung.ipynb)
* [20_Exceptions](http://nbviewer.jupyter.org/urls/gitlab.gwdg.de/jschulz1/math_prog/raw/master/lecture/20_Exceptions.ipynb)
* [21_GUI](http://nbviewer.jupyter.org/urls/gitlab.gwdg.de/jschulz1/math_prog/raw/master/lecture/21_GUI.ipynb)
* [22_performance](http://nbviewer.jupyter.org/urls/gitlab.gwdg.de/jschulz1/math_prog/raw/master/lecture/22_performance.ipynb)
* [23_Parallelisierung](http://nbviewer.jupyter.org/urls/gitlab.gwdg.de/jschulz1/math_prog/raw/master/lecture/23_Parallelisierung.ipynb)
* [24_git](http://nbviewer.jupyter.org/urls/gitlab.gwdg.de/jschulz1/math_prog/raw/master/lecture/24_git.ipynb)
......
%% Cell type:markdown id: tags:
# Aufgabe 1
1. Laden Sie die Datei `census.csv` die U.S. Population von 1790 bis 1990 in ihren Speicher und stellen Sie die Zahlen grafisch dar.
2. Interpolieren Sie die Kurve mit `interp1d` mit verschiedenen Methoden und stellen Sie die Ergebnisse grafisch dar.
3. Extrapolieren Sie mit Hilfe von kubischen Splines (dafür braucht man eine andere scipy-Funktion als `interp1d`) die Bevölkerungszahl im Jahr 2050. Da die Funktion ähnlich einer quadratischen Funktion ist macht eine Extrapolation mit kubischen/qudaratischen Funktionen *halbwegs* Sinn.
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
# Aufgabe 2
Für ein Polynom
$$ P(x) = x^n + p_1 x^{n-1} + \ldots + p_{n-1} x + p_n, $$
mit Leitkoeffizient 1 ist die Begleitmatrix $C_P$ definiert als
$$ \begin{pmatrix}
0 & 0 & 0 & \dots & 0 & -p_n \\
1 & 0 & 0 & \dots & 0 & -p_{n-1} \\
0 & 1 & 0 & \dots & 0 & -p_{n-2} \\
0 & 0 & 1 & \dots & 0 & -p_{n-3} \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & 0 & \dots & 1 & -p_1
\end{pmatrix}. $$
Ihr charakteristisches Polynom ist gerade $P$, ihre Eigenwerte sind daher die Nullstellen von $P$.
Schreiben Sie eine Funktion, die ein Polynom als als Argument hat und die Nullstellen zurückgibt. Stellen Sie dann das Polynom
$$ P(x)= x^5 - 4x^4 - 10x^3 + 40x^2 + 9x - 36$$
mit seinen Nullstellen zusammen graphisch dar.
*Anmerkung*: Man kann das Polynom sowohl als `poly1d`-Objekt als auch als SymPy-Polynom darstellen. Als Übung kann man auch beides implementieren. Die Berechnung der Eigenwerte sollte aber in beiden fällen numerisch erfolgen.
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
# Aufgabe 3
- Interpolieren Sie mit Polynomen die durch
$$f(x):=x^2\exp(-|x|)$$
gegebene Funktion an den Stellen `x = linspace(-5, 5, 13)`.
- Berechnen Sie approximativ den maximalen absoluten Fehler zwischen $f$ und ihrer Interpolierenden auf $[-5,5]$.
- Ändern Sie den Vektor der Stützstellen so, dass
$$ x_k = 5 \cos\left(\frac{2k+1}{26}\pi \right), \quad k=0, \dots , 12. $$
Berechnen Sie erneut den maximalen Fehler.
- Betrachten Sie auch die Stützstellen
$$ x_k = 5 \cos\left(\frac{2k+1}{100}\pi \right), \quad k=0, \dots , 49. $$
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
# Aufgabe 4
- Generieren Sie ein äquidistantes Gitter im Quadrat $[0,\, 1] \times [0,\, 1]$ mit 8 Punkten pro Achse und berechnen Sie die Funktion
$$ f(x,y) = \sin(4 \pi x) \cos(4 \pi y) $$
auf den Gitterpunken.
- Plotten Sie zuerst nur die Punkte.
- Interpolieren sie die Punkte auf einem feineres Gitter und erstellen Sie Grafiken mit `imshow` und `contour`.
- Beschriften Sie die Konturlinien von `contour`.
- Interpolieren sie die Punkte auf einem feineres Gitter und erstellen Sie Grafiken mit `Image` und `Contour`.
- Untersuchen Sie den Einfluss der verschiedenen Interpolationsmethoden.
%% Cell type:code id: tags:
``` python
```
......
%% Cell type:markdown id: tags:
# Aufgabe 1
1. Laden Sie die Datei `census.csv` die U.S. Population von 1790 bis 1990 in ihren Speicher und stellen Sie die Zahlen grafisch dar.
2. Interpolieren Sie die Kurve mit `interp1d` mit verschiedenen Methoden und stellen Sie die Ergebnisse grafisch dar.
3. Extrapolieren Sie mit Hilfe von kubischen Splines (dafür braucht man eine andere scipy-Funktion als `interp1d`) die Bevölkerungszahl im Jahr 2050. Da die Funktion ähnlich einer quadratischen Funktion ist macht eine Extrapolation mit kubischen/qudaratischen Funktionen *halbwegs* Sinn.
%% Cell type:code id: tags:
``` python
from numpy import loadtxt
import matplotlib.pyplot as plt
import holoviews as hv
hv.extension('bokeh')
date, pop = loadtxt('census.csv', delimiter=',', unpack=True)
plt.figure()
plt.plot(date, pop)
plt.show()
date, pop = loadtxt('../../data/census.csv', delimiter=',', unpack=True)
hv.Curve((date,pop), kdims='date', vdims='population')
```
%% Cell type:code id: tags:
``` python
from scipy.interpolate import interp1d
from numpy import linspace
plt.plot(date, pop, label='data')
hv.Curve((date, pop), label='data')
# Interpolations-Punkte
idate = linspace(min(date), max(date), 40)
pl = []
methods = ['nearest','linear','cubic']
for k, method in enumerate(methods):
f = interp1d(date, pop, kind=method, bounds_error=False)
ipop = f(idate)
plt.plot(idate, ipop, label=method)
pl.append(hv.Curve((idate, ipop), label=method))
plt.legend(loc=4)
plt.show()
hv.Overlay(pl)
```
%% Cell type:code id: tags:
``` python
from scipy.interpolate import UnivariateSpline
# Alle Interpolationen mit interp1d interpolieren aucht tatsaechlich nur,
# Extrapolation ist nicht möglich. Hier müssen wir uns für die Abschätzung
# etwas anderes einfallen lassen. Nutze UnivariateSpline:
f = UnivariateSpline(date, pop)
idate = linspace(min(date), 2050, 50)
plt.plot(idate, f(idate))
plt.show()
print(f(2050))
hv.Curve((idate, f(idate)))
```
%% Cell type:markdown id: tags:
# Aufgabe 2
Für ein Polynom
$$ P(x) = x^n + p_1 x^{n-1} + \ldots + p_{n-1} x + p_n, $$
mit Leitkoeffizient 1 ist die Begleitmatrix $C_P$ definiert als
$$ \begin{pmatrix}
0 & 0 & 0 & \dots & 0 & -p_n \\
1 & 0 & 0 & \dots & 0 & -p_{n-1} \\
0 & 1 & 0 & \dots & 0 & -p_{n-2} \\
0 & 0 & 1 & \dots & 0 & -p_{n-3} \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & 0 & \dots & 1 & -p_1
\end{pmatrix}. $$
Ihr charakteristisches Polynom ist gerade $P$, ihre Eigenwerte sind daher die Nullstellen von $P$.
Schreiben Sie eine Funktion, die ein Polynom als als Argument hat und die Nullstellen zurückgibt. Stellen Sie dann das Polynom
$$ P(x)= x^5 - 4x^4 - 10x^3 + 40x^2 + 9x - 36$$
mit seinen Nullstellen zusammen graphisch dar.
*Anmerkung*: Man kann das Polynom sowohl als `poly1d`-Objekt als auch als SymPy-Polynom darstellen. Als Übung kann man auch beides implementieren. Die Berechnung der Eigenwerte sollte aber in beiden fällen numerisch erfolgen.
%% Cell type:code id: tags:
``` python
from numpy import eye
from scipy.linalg import eigvals
def companion(p):
C = eye(p.order, k=-1)
C[:, -1] = -p.coeffs[:0:-1]
return C
def zeros(p):
return eigvals(companion(p))
```
%% Cell type:code id: tags:
``` python
from numpy import poly1d, linspace
import matplotlib.pyplot as plt
p = poly1d([1, -4, -10, 40, 9, -36])
z = zeros(p)
print(p)
# `zeros` liefert komplexe Zahlen, was für das Plotten problematisch ist.
# Tatsächlich sind die Nullstellen hier aber reell. Konsistenzcheck:
print('Maximaler Imaginärteil = ', max(abs(z.imag)))
z = z.real
x = linspace(min(z) - 0.2, max(z) + 0.2, 100)
plt.plot(x, p(x))
plt.plot(z, p(z), 'o')
plt.show()
hv.Curve((x,p(x))) * hv.Points((z, p(z))).options(size=10)
```
%% Cell type:code id: tags:
``` python
# Alternativ: SymPy-Lösung
from sympy import poly, symbols
from numpy import array, eye
from numpy.linalg import eigvals
def companion(p):
C = eye(p.degree(), k=-1)
C[:, -1] = -array(p.coeffs())[:0:-1]
return C
def zeros(p):
return eigvals(companion(p))
x = symbols('x')
expr = x**5 - 4*x**4 - 10*x**3 + 40* x**2 + 9*x -36
p = poly(expr)
z = zeros(p)
z = z.real
```
%% Cell type:code id: tags:
``` python
# Zum Auswerten des Polynoms wird hier Theano verwendet
from sympy.printing.theanocode import theano_function
# Zum Auswerten des Polynoms wird hier lambdify verwendet
from sympy import lambdify
f = theano_function([x], [p.as_expr()], dims={x: 1})
f = lambdify([x], p.as_expr())
xs = linspace(min(z) - 0.2, max(z) + 0.2, 100)
plt.plot(xs, f(xs))
plt.plot(z, f(z), 'o')
plt.show()
hv.Curve((xs, f(xs))) * hv.Points((z, f(z))).options(size=10)
```
%% Cell type:code id: tags:
``` python
p.as_expr()
```
%% Cell type:markdown id: tags:
# Aufgabe 3
- Interpolieren Sie mit Polynomen die durch
$$f(x):=x^2\exp(-|x|)$$
gegebene Funktion an den Stellen `x = linspace(-5, 5, 13)`.
- Berechnen Sie approximativ den maximalen absoluten Fehler zwischen $f$ und ihrer Interpolierenden auf $[-5,5]$.
- Ändern Sie den Vektor der Stützstellen so, dass
$$ x_k = 5 \cos\left(\frac{2k+1}{26}\pi \right), \quad k=0, \dots , 12. $$
Berechnen Sie erneut den maximalen Fehler.
- Betrachten Sie auch die Stützstellen
$$ x_k = 5 \cos\left(\frac{2k+1}{100}\pi \right), \quad k=0, \dots , 49. $$
%% Cell type:code id: tags:
``` python
from numpy import exp
def f(x):
return x**2 * exp(-abs(x))
```
%% Cell type:code id: tags:
``` python
from scipy.interpolate import BarycentricInterpolator
import matplotlib.pyplot as plt
def interpolation_and_error(f, x):
p = BarycentricInterpolator(x, f(x))
# Zur Berechnung des Fehlers werden beide Funktionen auf einem
# feineren Gitter ausgewertet.
z = linspace(x[0], x[-1], 200)
return p, max(abs(f(z) - p(z)))
def plot_with_interpolation(f, p, x):
def plot_with_interpolation(f, p, x, title):
plot_x = linspace(x[0], x[-1], 200)
plt.plot(plot_x, f(plot_x), 'b--', linewidth=3, label='Funktion')
plt.plot(plot_x, p(plot_x), 'r', linewidth=3, label='Interpolierende')
plt.plot(x, f(x), 'o', label='Stützstellen')
plt.legend(loc='best')
plt.grid('on')
return (hv.Curve((plot_x, f(plot_x)), label='Funktion') * hv.Curve((plot_x, p(plot_x)), label='Interpolierende') * hv.Points((x, f(x)), label='Stützstellen')).relabel(title)
def make_plot(x):
p, err = interpolation_and_error(f, x)
plot_with_interpolation(f, p, x)
plt.title('err = {}'.format(err))
plt.show()
return plot_with_interpolation(f, p, x, 'err = {}'.format(err))
```
%% Cell type:code id: tags:
``` python
from numpy import linspace
make_plot(linspace(-5, 5, 13))
```
%% Cell type:code id: tags:
``` python
from numpy import ogrid, cos, pi, arange
def sample_points(N):
return 5 * cos(pi * (2*arange(N)+1) / (2*N))
```
%% Cell type:code id: tags:
``` python
make_plot(sample_points(13))
```
%% Cell type:code id: tags:
``` python
make_plot(sample_points(50))
```
%% Cell type:markdown id: tags:
# Aufgabe 4
- Generieren Sie ein äquidistantes Gitter im Quadrat $[0,\, 1] \times [0,\, 1]$ mit 8 Punkten pro Achse und berechnen Sie die Funktion
$$ f(x,y) = \sin(4 \pi x) \cos(4 \pi y) $$
auf den Gitterpunken.
- Plotten Sie zuerst nur die Punkte.
- Interpolieren sie die Punkte auf einem feineres Gitter und erstellen Sie Grafiken mit `imshow` und `contour`.
- Beschriften Sie die Konturlinien von `contour`.
- Interpolieren sie die Punkte auf einem feineres Gitter und erstellen Sie Grafiken mit `Image` und `Contour`.
- Untersuchen Sie den Einfluss der verschiedenen Interpolationsmethoden.
%% Cell type:code id: tags:
``` python
from numpy import sin, cos, pi, array, ogrid
from numpy import sin, cos, pi, array, mgrid
from numpy.random import random
def f(x, y):
return cos(4*pi*x) * sin(4*pi*y)
X, Y = ogrid[0:1:8j, 0:1:8j]
X, Y = mgrid[0:1:8j, 0:1:8j]
Z = f(X, Y)
```
%% Cell type:code id: tags:
``` python
import matplotlib.pyplot as plt
from numpy import broadcast_arrays
plt.scatter(*broadcast_arrays(X, Y), c=Z)
plt.axis('scaled')
plt.xlim((0, 1))
plt.ylim((0, 1))
plt.figure()
plt.imshow(Z.T, interpolation='none', extent=(0, 1, 0, 1))
plt.show()
hv.Scatter((X, Y,Z), vdims=['y', 'scaled']).options(color_index=2, size=10) + hv.Image(Z, bounds=(0, 0, 1, 1))
```
%% Cell type:code id: tags:
``` python
from scipy.interpolate import interp2d
from numpy import linspace
def make_plots(kind):
p = interp2d(X, Y, Z, kind=kind)
# Aufrufen von interp2d-Objekten ist ungewöhnlich insofern, als dass für beide Argumente 1d-Arrays
# erwartet werden und das Broadcasting implizit gemacht wird.
plot_axis = linspace(0, 1, 100)
vals = p(plot_axis, plot_axis)
plt.figure()
plt.imshow(vals.T, interpolation='none', extent=(0, 1, 0, 1))
img = hv.Image(vals, bounds=(0, 0, 1, 1))
plt.figure()
plt.contourf(vals)
c = plt.contour(vals, color='black')
plt.clabel(c, inline=1)
```
%% Cell type:code id: tags:
``` python
make_plots('linear')
plt.show()
return img * hv.operation.contours(img).options(show_legend=False)
```
%% Cell type:code id: tags:
``` python
make_plots('cubic')
plt.show()
make_plots('linear')+make_plots('cubic')+make_plots('quintic')
```
%% Cell type:code id: tags:
``` python
make_plots('quintic')
plt.show()
```
......