12.6 Berechnung von Reihen und Konvergenzverhalten

Ziel:

Ziel der Übung ist die Berechnung von endlichen und unendlichen Reihen in MATLAB.

Anzufertigen sind die MATLAB-Funktionen reihe und eulernum. Die Funktion reihe löst die Reihenentwicklung einiger Funktionen, und die Funktion eulernum berechnet die Eulerschen Zahlen, die für eine spezielle Reihe notwendig sind.

Darüber hinaus benötigt man ein MATLAB-Skript uebreihe, welches die Funktionen testet und die Ergebnisse in einfachen Plots darstellt. Exemplarische Plots sind nach der Beschreibung der Übung eingefügt.

Voraussetzung:

Die Voraussetzung für die Übung sind Kenntnisse über die grundlegenden Befehle in MATLAB und über die Erstellung von Funktionen und deren Verwendung. Außerdem braucht man einfache Plotbefehle, wie plot, legend, xlabel, ylabel, title.

Als einfaches Beispiel für eine unendliche Reihe wird hier die Reihendarstellung der Funktion $\sin(x)$ verwendet. Diese ist gegeben durch

$$\sin(x) = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} + \dots$$ (12.7)

$$= \sum_{k=0}^{\infty} (-1)^k \frac{x^{2k+1}}{(2k+1)!} \; ,$$ (12.8)

welche für alle $\vert x\vert<\infty$ konvergiert.

In allen numerischen Programmiersprachen kann man solche unendlichen Reihen natürlich nur näherungsweise auswerten, da die Summation von unendlich vielen Termen unendlich lange dauern würde. Daher gibt man sich mit der endlichen Teilsumme

$$\sin(x) \approx \sum_{k=0}^{n} (-1)^k \frac{x^{2k+1}}{(2k+1)!} \; ,$$ (12.9)

zufrieden. Abhängig von der gewählten Reihe und von den Werten für $x$, konvergieren sie langsamer, schneller oder eben gar nicht zum gewünschten Wert.

Summen dieser Art kann man in MATLAB sehr gut ohne die Verwendung von for-Schleifen programmieren. Nützlich sind dabei folgende Informationen:

• Die Berechnung von $k!$ (Fakultät) kann mit Hilfe der $\Gamma$-Funktion durchgeführt werden

  $$\Gamma\left( n+1 \right) = n! \qquad \left( n = 0,1,2,\dots\right) \; ,$$ (12.10)

  
  
  wobei dafür in MATLAB die Funktion gamma(k) zur Verfügung steht. Diese funktioniert natürlich wieder für Arrays von $k$-Werten. Für große $k > 171$ liefert $\Gamma(k)$ den Wert unendlich, da das Ergebnis jenseits der Darstellungsgrenze für den Datentyp double liegt.
• Äußerst praktisch für die Berechnung solcher Reihen ist der Befehl meshgrid. Damit kann man aus einem Vektor ${\bf x}$ mit allen $x$-Werten und einem Vektor ${\bf k}$ mit allen $k$-Werten zwei gleich große Matrizen

   x = [-2:2]; k = [0:5];
   [xx,kk] = meshgrid(x,k);
  

  erzeugen

  $$xx = \begin{bmatrix}-2 & -1 & 0 & 1 & 2 -2 & -1 & 0 & 1 & 2 \\... ... 3 & 3 & 3 & 3 4 & 4 & 4 & 4 & 4 5 & 5 & 5 & 5 & 5 \end{bmatrix} \qquad .$$ (12.11)

  
  
  Mit dem Befehl A=xx.^(2*kk+1) kann man nun z.B. $x^{2k+1}$ für alle Werte von $x$ und $k$ gleichzeitig ausrechnen. Die Summe über alle $k$-Werte kann man dann mit Hilfe des sum(A,1)-Befehles erhalten.
• Sehr einfach kann die gesamte Reihenauswertung auch mit Hilfe der Matrizenmultiplikation durchgeführt werden. Aus der linearen Algebra ist bekannt, dass für einen $(1 \times n)$-Vektor ${\bf f}$ (Zeilenvektor) und eine $(n \times m)$-Matrix ${\bf A}$ gilt

  $$\{{\bf f}{\bf A}\}_l = \sum_{k} f_k A_{kl} \; .$$ (12.12)

  
  
  Erzeugt man nun für alle Werte von $k$ den Vektor ${\bf f}= (-1)^k / (2k+1)!$, kann man mit ${\bf f}{\bf A}$ (in MATLAB f*A) für alle Werte von $x$ gleichzeitig den jeweiligen Wert der Teilsumme erhalten.

Der Aufruf der MATLAB-Funktion reihe soll folgendermaßen aussehen

 [r,u,a] = reihe(typ,x,n)

wobei typ eine Schaltvariable (switch), x ein Array von $x$-Werten und n der Maximalwert für die $k$-Werte ist. Ausgegeben werden sollen die Teilsummen r, die Werte der Funktion u (hier z.B. $\sin(x)$) und eine eventuelle analytische Lösung a für die endlichen Teilsummen der Reihe.

Alle Outputwerte müssen die gleiche Größe wie $x$ haben. Dazu merkt man sich am Anfang des Programms die Größe von $x$ (size), macht aus dem Array $x$ einen Vektor (colon). Am Ende des Programms kann man dann mit reshape sicherstellen, dass die richtige Größe zurückgegeben wird.

1. Programmieren Sie entsprechend der Angabe die Reihe für die Sinusfunktion (typ='sin'). Für $a$ geben Sie lauter NaN zurück.

   Stellen Sie mit Hilfe des Skripts uebreihe einen Vergleich zwischen der Teilsumme der Reihe und der Funktion für verschiedene Werte von $n$ dar.

2. Die Teilsumme der Geometrischen Reihe ist gegeben durch

   $$r=\sum_{k=1}^{n} x^{k-1} = \left\{ \begin{array}{r@{\quad:\quad}l} (1 - x^n) / (1 - x) & x \neq 1 n & x = 1 \end{array} \right. \; .$$ (12.13)

   
   
   Die unendliche geometrische Reihe konvergiert für alle $\left\vert x \right\vert < 1$

   $$u=\sum_{k=1}^{\infty} x^{k-1} = \left\{ \begin{array}{r@{\quad:\q... ... x\vert < 1 \infty & x \geq 1 NaN & x \leq -1 \end{array} \right. \; .$$ (12.14)

   
   
   Berechnen Sie für den typ 'geom' den Wert r der geometrischen Reihe 12.13, in a den analytischen Wert aus 12.13 und in u den Wert der unendlichen Reihe 12.14. Die Fallunterscheidungen können Sie mit Hilfe der logischen Indizierung durchführen.

   Probieren Sie in uebreihe das Konvergenzverhalten aus und stellen Sie die Ergebnisse für r, a und u gemeinsam als Funktion von x dar.

3. Eine der Hyperbelfunktionen, ${\rm sech}(x)=1/{\rm cosh}(x)$, kann durch folgende asymptotische Reihe dargestellt werden

   $${\rm sech}(x) = 2 \sum_{k=0}^{\infty} (-1)^k e^{-(2k+1) \left\vert x\right\vert} \; ,$$ (12.15)

   
   
   die nun für große Werte von $x$ rasch konvergiert und um Null Schwierigkeiten macht . Als Typ verwenden Sie 'sechasymp'.

   Geben Sie in r die Werte der Teilsumme, in u die Werte der Funktion sech und in a einen Vektor mit NaN zurück.

4. Die Taylorreihenentwicklung für ${\rm sech}$ ist gegeben durch

   $${\rm sech}(x) = 1 - \frac{x^2}{2} + \frac{5x^4}{24} - \frac{61x^6}{720} + \dots$$ (12.16)

   $$= \sum_{k=0}^{\infty} \frac{{\rm E}_k}{(2k)!} x^{2k} \; ,$$ (12.17)

   
   
   welche für alle $\vert x\vert<\pi/2$ konvergiert (typ='sech'). Außerhalb des Konvergenzintervalls geben Sie NaN zurück.

   Die Eulerschen Zahlen ${\rm E}_k$ sind durch folgende Reihe definiert

   $${\rm E}_k = \frac{(-1)^k 2^{2k+2} (2k)!}{\pi^{2k+1}} \sum_{l=1}^{\infty} \frac{(-1)^{l-1}}{(2l-1)^{2k+1}} \; ,$$ (12.18)

   
   
   die im Programm

    e = eulernum(k,lmax)
   

   berechnet werden sollen. Hier wird nun über $l$ summiert und der ${\bf k}$-Vektor spielt die Rolle, die vorher $x$ gespielt hat, d.h. für jeden Wert von $k$ bekommt man einen Wert von ${\rm E}_k$. Als Defaultwert für lmax kann man 10 verwenden. Da alle Eulerzahlen ganzzahlig sind, empfiehlt es sich, am Ende der Funktion den Befehl round zu verwenden. Damit werden leichte numerische Ungenauigkeiten ausgeglichen, und man muss lmax nicht größer wählen. Zur Kontrolle hier die Werte von ${\rm E}_0$ bis ${\rm E}_6$, die sich als 1, -1, 5, -61, 1385, -50521, 2702765 ergeben. Diese sollten Sie beim Aufruf von eulernum([0:6]) erhalten.

   Vergleichen Sie im Skript uebreihe die Ergebnisse der beiden Reihenentwicklungen und stellen Sie sie graphisch dar.