12.14 Gaußfunktion

Ziel:

Bei der vorliegenden Übung wird nochmals die Erstellung von Funktionen geübt. Darüber hinaus wird hier auch Wert auf die graphische Darstellung in zwei und drei Dimensionen gelegt. Abzugeben sind die Unterprogramme
gauss.m, gauss1d.m, gauss2d.m und gausstest.m.

Voraussetzung:

Voraussetzung für die Übung ist das Kapitel über Programmeinheiten 7 und über Graphik 11. Als kleine Einführung zum Plotten können die Files htplot2d.m und htplot2da.m der verwendet werden.

In der Wahrscheinlichkeitsrechnung und auch in der Physik hat die Normalverteilung bzw. Gaußverteilung eine große Bedeutung. Sie ist definiert durch

$$g(x)= \frac{1}{\sqrt{2\pi}\sigma} \exp\left\{-\frac{(x-x_{0})^2}{2 \sigma^2}\right\} \; ,$$ (12.54)

wobei $x_{0}$ und $\sigma$ die Parameter der Verteilung sind. In $x_{0}$ liegt sowohl das Maximum als auch das Symmetriezentrum, und $\sigma$ ist der Abstand von diesem Zentrum zu den Wendepunkten. Wie für jede Wahrscheinlichkeitsverteilung gilt

$$\int_{-\infty}^{\infty} g(x) dx = 1 \; .$$ (12.55)

In der Physik verwendet man häufig auch eine Summation über mehrere Gaußfunktionen mit jeweils unterschiedlichen Parametern. Dies hat den Sinn, dass man gleichzeitig mehrere Maxima darstellen kann. Allgemein kann man die Funktion dann definieren als

$$g(x)= \sum_{k=1}^{n} \frac{f_k}{\sqrt{2\pi}\sigma_k} \exp\left\{-\frac{(x-x_{0k})^2}{2 \sigma_k^2}\right\} + g_0 \; ,$$ (12.56)

wobei $x_{0k}$ und $\sigma_k$ die gleiche Bedeutung wie vorher haben. Der Parameter $f_k$ ist ein Multiplikator, um Peaks unterschiedlicher Höhe darstellen zu können. Die Summation erfolgt über alle $n$ Werte der Parameter $x_{0k}$, $\sigma_k$ und $f_k$. Zusätzlich kann die ganze Funktion noch um den Wert $g_0$ vertikal verschoben sein. Für $g_0=0$ gilt natürlich

$$\int_{-\infty}^{\infty} g(x) dx = \sum_{k=1}^{n} f_k \; .$$ (12.57)

In zwei Dimensionen kann die Gaußverteilung folgendermaßen definiert werden

$$g(x,y) = \frac{1}{2\pi \sigma_{x}\sigma_{y}} \exp\left\{ -\frac{(x-x_{0})^2}{2 \sigma_{x}^2} - \frac{(y-y_{0})^2}{2 \sigma_{y}^2} \right\} \; ,$$ (12.58)

wobei nun $x_{0}$, $y_{0}$, $\sigma_x$ und $\sigma_y$ die Parameter der Verteilung sind. Auch hier gilt die Normierungsbedingung

$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} g(x,y) dx dy = 1 \; .$$ (12.59)

Die Summation über mehrere Peaks kann natürlich auch hier erfolgen,

$$g(x,y)= \sum_{k=1}^{n} \frac{f_k}{2\pi \sigma_{xk} \sigma_{yk} } ... ...\sigma_{xk}^2} - \frac{(y-y_{0k})^2}{2 \sigma_{yk}^2} \right\} + g_0 \; ,$$ (12.60)

wobei die Parameter die analoge Bedeutung wie in 12.56 haben. Für $g_0=0$ gilt natürlich

$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} g(x,y) dx dy = \sum_{k=1}^{n} f_k \; .$$ (12.61)

1. Schreiben Sie nun eine Funktion gauss1d.m, welche die Formel 12.56 für beliebige numerische Arrays von $x$ auswertet. Die Ausgabegröße sollte ein Feld $g$ sein, das die gleiche Größe und Dimension wie $x$ hat. Die Parameter $x_{0k}$, $\sigma_k$ und $f_k$ können ebenfalls in numerischen Arrays gespeichert sein, welche die gleiche Größe haben müssen. Der Parameter $g_0$ muss natürlich ein Skalar sein.

   • Die Deklaration der Funktion sollte so aussehen.
     function [g] = gauss1d(x,x0,s,f,g0)
     x0 enthält alle $x_{0k}$-Werte, s alle $\sigma_{k}$-Werte, und f alle $f_{k}$-Werte.
   • Setzen Sie die Defaultwerte $x_{0}=[-1,1]$, $\sigma=[0.5,1]$, $f=[1,0.5]$ und $g_0=0$ (nargin).
   • Überprüfen Sie die Eingabegrößen.
   • Initialisieren Sie $g$ als Feld mit lauter Nullen in der Größe von $x$.
   • Bilden Sie die Summe über alle Werte ( $k = 1\dots n$) der Parameter mit Hilfe einer einzigen for-Schleife. Vergessen Sie dabei nicht, dass $x$ ein Array sein kann.
   • Probieren Sie die Funktion mit Hilfe des Skripts gausstest.m aus und machen Sie einen Plot
     plot(x,g)
     Wählen Sie den $x$-Vektor in einem vernünftigen Bereich, damit die Funktion $g(x)$ sinnvoll dargestellt wird (linspace). Beschriften Sie den Plot mit Achsenbeschriftung und Titel.

2. Schreiben Sie nun eine Funktion gauss2d.m, welche die Formel 12.60 auswertet. Die Inputgrößen $x$ and $y$ können entweder gleich große Arrays sein oder eine der beiden ist ein Array und die andere ein Skalar. Im ersten Fall (gleich groß) wird die Funktion für alle Punkte $(x_i,y_i)$ berechnet, im zweiten Fall (ein Skalar) für alle Punkte $(x_i,y_1)$ oder $(x_1,y_i)$. Die Ausgabegröße sollte ein Feld $g$ sein, das die gleiche Größe und Dimension wie das Array $x$ oder $y$ hat. (Aufpassen, da $x$ oder $y$ auch ein Skalar sein können!) Die Parameter $x_{0k}$, $\sigma_{xk}$, $y_{0k}$, $\sigma_{yk}$ und $f_k$ können ebenfalls in numerischen Arrays gespeichert sein, welche die gleiche Größe haben müssen. Der Parameter $g_0$ muss natürlich wieder ein Skalar sein.

   • Die Deklaration der Funktion sollte so aussehen.
     function [g] = gauss2d(x,y,x0,sx,y0,sy,f,g0)
   • Setzen Sie die Defaultwerte für zwei Peaks Ihrer Wahl.
   • Überprüfen Sie die Eingabegrößen.
   • Initialisieren Sie $g$ als Feld mit lauter Nullen in der Größe von $x$ oder $y$.
   • Bilden Sie nun die Summe in 12.60 mit einer for-Schleife über alle Werte $k = 1\dots n$.
   • Probieren Sie die Funktion mit Hilfe des Skripts gausstest.m aus und machen Sie eine dreidimensionale Darstellung
     surf(xx,yy,g) oder
     mesh(xx,yy,g).
     Wählen Sie den $x$-Vektor und den $y$-Vektor in einem vernünftigen Bereich und bilden Sie mit meshgrid Arrays für die Darstellung der Funktion.
   • Beschriften Sie den Plot und probieren Sie sinnvolle Darstellungen in Subplots einer Figure. Zum Beispiel sollte man auch Konturplots ausprobieren.

3. Eine numerische Integration kann in MATLAB mit Hilfe der Funktion quadl durchgeführt werden. Ein Aufruf zum Auswerten des 1-D Integrals kann so aussehen,
   flaeche = quadl('gauss',x1,x2,tol,trace,parameter),
   wobei hier gauss die Funktion ist, über die von x1 bis x2 integriert wird. Für die Genauigkeit tol und für den Wert trace können zum Testen leere Felder [] eingegeben werden. Danach folgen die Parameter der Gaußfunktion x0, s, f und y0.
   • Schreiben Sie am Ende von gausstest.m den Code für die Integration und geben Sie das Ergebnis formatiert aus. Überprüfen Sie es mit der analytischen Lösung 12.57. Was passiert, wenn $g_0$ von Null verschieden ist?
   • Probieren Sie es mit verschiedenen Parametern aus und überprüfen Sie 12.57. Bedenken Sie dabei, dass die Integrationsgrenzen weit genug entfernt von den Maxima sind, damit 12.57 näherungsweise erfüllt sein kann. Natürlich dürfen sie aber nicht $\pm \infty$ sein.
   • Machen Sie das Gleiche für die 2-dimensionale Gaußfunktion.

4. • Schreiben Sie die Funktion gauss.m, die auf zwei Arten aufgerufen werden kann,
     [g] = gauss(a,x) oder
     [g] = gauss(a,x,y).
   • Aus ihr werden je nach Anzahl der Inputparameter (2 oder 3) die 1-D oder 2-D Funktionen gauss1d oder gauss2d aufgerufen. Der Vektor a soll dabei alle Parameter enthalten, also
     $[x_{0k},\sigma_{k},f_{k},g_{0}]$ in 2-D, oder
     $[x_{0k},\sigma_{xk},y_{0k},\sigma_{yk},f_{k},g_{0}]$ in 3-D.
   • Diese Vektoren haben die Länge $3n+1$ bzw. $5n+1$, wobei $n$ die Anzahl der Peaks ist. Vor dem Aufruf von gauss1d oder gauss2d müssen Sie auf die dort benötigten Inputparameter aufgespalten werden. Bestimmen Sie daher zuerst die Länge von a, daraus dann n, und damit sollte die Aufspaltung klar sein. Sollte die Länge von a nicht stimmen, soll eine error-Meldung ausgegeben werden.
   • Testen Sie die neue Routine im Skript gausstest.m mit zwei möglichst sinnvollen Aufrufen.

Geben Sie gauss.m, gauss1d.m, gauss2d.m und das Skript gausstest.m ab. Verwenden Sie dabei wieder !uebungsabgabe.