Programmierpraktikum im WS 2011/2012

Denny Otten, V5-142, dotten@math.uni-bielefeld.de

Pflichtaufgaben 2: Fraktale mit Matlab

Dieser Zettel besteht aus drei Teilen, es müssen zur erfolgreichen Teilnahme am Kurs alle Pflichtaufgaben aus mindestens zwei der drei Teile erfolgreich bearbeitet werden.

1 Lindenmayer-Systeme – Pflichtaufgaben Teil 1

Als erstes betrachten wir einfache Lindenmayer-Systeme zur Erzeugung von fraktaler Kurven. Dazu sei Ω ein Alphabet von Zeichen, V ⊂ Ω die Menge aller Variablen und C ⊂ Ω die Menge aller Konstanten, Ω = V ˙∪C. Ferner sei W die Menge aller Worte über Ω und P eine Abbildung P : V → W die Produktionsregeln. Mit diesen Bezeichnungen ist ein Lindenmayer-System gegeben als ein 4-Tupel G = {V, C, ω, P }, wobei ω ein Wort aus W ist. Die Regeln werden iterativ auf jede Variable aus ω angewandt, wobei alle Produktionsregeln, gleichzeitig angewandt werden.

Beispiel. Fibonacci-Zahlen: Ω = {A, B}, V = {A, B}, C = {}, P (A) = B, P (B) = AB, ω = A, damit:

\begin{matrix}n = 0: & A & n = 1: & B & n = 2: & AB \\ n = 3: & BAB & n = 4: & ABBAB & n = 5: & BABABBAB \end{matrix}

Man kann Lindenmayer-Systeme mit Hilfe sogenannter Turtle-Grafik interpretieren. Eine abgespeckte, für unsere Zwecke völlig ausreichende Version liefert die folgende Tabelle:

Symbol Aktion
F Vorwärtsbewegung mit Zeichenstift
+ Linksdrehung um vorgegebenem Winkel δ
− Rechtsdrehung um vorgegebenem Winkel δ

Implementierung in Matlab:

[figure]

  function [xdata,ydata]=tgr(txt,x0,y0,phi0,delta)
    N = length(find(txt=='F'));
    xdata = zeros(1,N); ydata = zeros(1,N);
    phi = phi0;
    xdata(1) = x0; ydata(1) = y0;
    n = 1;
    for j = txt
       switch j
       case '+'
          phi = phi+delta;
       case '-'
          phi = phi-delta;
       case 'F'
          xdata(n+1) = xdata(n)+cos(phi);
          ydata(n+1) = ydata(n)+sin(phi);
          n = n+1;

       otherwise

          error('Rule not specified.');

       end
    end
  end

[figure]

Beispiel (Koch Kurve). Die Koch Kurve hast Du schon auf dem dritten Blatt kennengelernt. Man kann sie auch mit dem Lindenmayer-System (V, C, ω, P), gegeben durch

V = {F }, C = {+, −}, P (F) = F + F − −F + F, ω = F, Ω = {F, +, −}

erzeugen. Dafür setzt man den Winkel δ auf \frac{π}{3} , und erhält verschiedene Iterationsstufen der Kochkurve:

[figure]

  function txt=lm(omega,produ,depth)
    txt = omega;
    pr = {'F',produ;'+','+';'-','-'};
    for j = 1:depth
       txtn = [];
       for t = txt
          txtn = [txtn,pr{strcmp(pr(:,1),t),2}];

       end

       txt = txtn;

    end
  end

[figure]

Dieses führt man dann wie folgt aus (dabei soll die Iterationstiefe 4 sein):

>> txt=lm('F','F+F--F+F',4);
>> [x,y]=tgr(txt,0,0,0,pi/3);
>> plot(x,y), axis equal;

Wie muss ω gewählt werden um die Koch-Schneeflocke zu bekommen?

Aufgabe 1 (Pflichtaufgabe). Versuche die Programme tgr und lm zu verstehen und kommentiere den Programmcode. Zeichne die Minkowski-Kurve, die durch folgendes Lindenmayer-System gegeben ist: Ω = \{F, +, −\}, V = \{F \}, C = \{+, −\}, P (F) = F + F − F − F F + F + F − F, ω = F, der Winkel soll δ = \frac{π}{2} sein.

Aufgabe 2 (Pflichtaufgabe). Ändere die Programme so ab, dass Du auch die Drachenkurve, die durch das Lindenmayer-System mit Ω = {F, G, +, −}, V = {F, G}, C = {+, −}, P (F) = F + G+ , P (G) = −F − G, ω = F gegeben ist, erzeugen kannst. Der Winkel ist dabei δ = \frac{π}{2} und beim Zeichnen in der Turtle-Grafik soll G genau wie F behandelt werden! Erzeuge die Drachenkurve bis zur 12ten Iterationstiefe.

2 Juliamengen und Apfelmännchen – Pflichtaufgaben Teil 2

Fraktale lassen sich nicht nur –wie oben gesehen– geometrisch erzeugen, sondern auch duch Iteration von Funktionen (d.h. wiederholte Anwendung einer Funktion). Dazu startet man mit einem Punkt z_{0} und führt die folgende Vorschrift aus

z_{n+1} = f(z_{n}), n = 0, 1, 2, …

Dabei kann man sich n als diskrete Zeit (etwa Jahre) vorstellen. Besonders interessant ist es, wenn die Funktion f von ℂ → ℂ abbildet. Schon die quadratische Funktion

q_{c} : ℂ → ℂ, z ↦ z^{2} + c, c ∈ ℂ ein Parameter

zeigt ein sehr überraschendes Verhalten. Durch Iteration mit dieser Funktion wird das berühmte Apfelmännchen-Fraktal erzeugt. Wir werden im Folgenden sehen, auf welche Weise das geschieht.

Aufgabe 3. • Wende die Abbildung q_{c} für c = 0, c = −1.25, und c = 0.7i jeweils mehrfach auf die Punkte z_{0} = −0.1 + i, z_{0} = −0.3 − 0.5i, z_{0} = −0.7i an und stelle die Folge der Bilder {z_{n}}_{n∈ℕ} mit der Funktion plot als Punkte in der komplexen Ebene dar. Verbinde die Punkte zusätzlich mit Linien um den Weg des Anfangspunktes unter der Abbildung q_{c} gut ablesen zu können.

• Schreibe eine reelle Version qr der Funktion q_{c}, d.h. definiere eine Funktion qr_{cr,ci} : ℝ^{2} → ℝ^{2}, die als Eingabe den Realteil und den Imaginärteil des Anfangswertes erhält und statt mit dem Parameter c mit dem Realteil cr = Re(c) und dem Imaginärteil ci = Im(c) arbeitet.

Julia und Gefangenenmenge:

Alle Startpunkte z_{0}, für die die Folge {z_{n}}_{n∈ℕ} mit

z_{n+1} = q_{c}(z_{n}) = z_{n}^{2} + c, n = 0, 1, 2, …

für alle Zeiten beschränkt bleibt,bilden die Gefangenenmenge G_{c} := {z_{0} : lim sup_{n→∞} |z_{n}| < ∞}. Der Rand dieser Menge J_{c} := ∂G_{c} wird nach dem französischen Mathematiker Gaston Julia (1893–1978) Juliamenge genannt. Das Aussehen der Juliamenge bzw. der Gefangenenmenge ist je nach Wahl des Parameters c ∈ ℂ sehr verschieden.

Aufgabe 4. • Definiere eine Funktion flucht, die für eine Zahl z_{0} = x_{0} + iy_{0} ∈ ℂ die Anzahl n der Iterationen ausgibt, nach denen für ein Folgenglied z_{n} gilt |z_{n}| > 2, (Man kann zeigen, dass die Folge {z_{n}}_{n∈ℕ} für fast alle Punkte z_{0} mit |z_{0}| > 2 divergiert).

• Eingabeparameter sollten der Realteil x_{0} und Imaginärteil y_{0} sein. Da man nicht unendlich viele Iterationsschritte testen kann, sollte die Iteration nach einer bestimmten Anzahl von Schritten abgebrochen werden (z.B. 25 Schritte). Die so definierte Funktion bildet also ℝ × ℝ auf ℕ ab.

• Teste diese Funktion mit den in Aufgabe 3 gegebenen Startwerten z_{0} und Parametern c.

Das Bild der Gefangenenmenge erhält man, wenn man nicht nur einen einzigen Anfangswert z_{0} ∈ ℂ betrachtet, sondern für jeden solchen Anfangswert z_{0} = x_{0} + iy_{0} aus einen bestimmten Gebiet das Ergebnis der Funktion flucht farblich kodiert zeichnet.

Aufgabe 5. Benutze die Funktion pcolor, imagesc oder surf ( pcolor(x,y, C), shading flat,
imagesc(x,y,C) oder surf(x,y, C), shading flat) um die Funktion flucht in einem geeigneten Be-

reich x_{0} = x_{min}, ..., x_{max}, y_{0} = y_{min}, ..., y_{max} in ℝ × ℝ zu zeichnen. Probiere dazu verschiedene Werte des Parameters c aus, etwa

−1, 0.3603 + 0.1i, −0.122 + 0.745i, i, 0.113 − 0.6704i.

Welches ist die Gefangenenmenge G_{c}?

Aufgabe 6 (Pflichtaufgabe). Schreibe eine Funktion flucht reell, die das gleiche wie die Funktion flucht in reeller Arithmetik tut. Werte dazu die Operationen wie Multiplikation, Addition, etc. in Realund Imaginärteil getrennt aus. Zeichne wie in Aufgabe 5. (Zusatz: Vektorisiere die Funktionen flucht und flucht reell und vergleiche die Laufzeiten.)

Mandelbrotmenge:

Das Apfelmännchen bzw. die Mandelbrotmenge, benannt nach dem Mathematiker Benoît Mandelbrot (1924–), besteht aus den Parameterwerten c ∈ ℂ, für die der Startwert z_{0} = 0 im Laufe der Iteration beschränkt bleibt. Man kann zeigen, dass die Mandelbrotmenge Teilmenge des Quadrates {c ∈ ℂ : |Re(c)| ≤ 2, |Im(c)| ≤ 2} enthalten ist. Das Apfelmännchen entsteht also ähnlich wie die Juliamenge: Um zu bestimmen, ob ein Punkt z ∈ ℂ in der Mandelbrotmenge liegt, startet man die Iteration z_{n+1} = z_{n}^{2}+c bei z_{0} = 0 und rechnet solange, bis |z_{n}| > 2 ist. Jeder Punkt c ∈ ℂ wird nach der Anzahl der Schritte eingefärbt, nach denen er den Kreis verlässt. Wenn der Punkt nach einer grossen Anzahl von Schritten (z.B. 25) den Kreis immer noch nicht verlassen hat, wird die Iteration abgebrochen.

Aufgabe 7 (Pflichtaufgabe). Schreibe wie oben Funktionen mandelbrot und mandelbrot reell, die für eine Zahl c ∈ ℂ die Anzahl der Iterationen n ausgibt, nach denen |z_{n}| > 2 gilt. Benutze beide Funktionen um die Mandelbrotmenge (d.h. das Apfelmännchen) zu zeichnen. (Zusatz: Vektorisiere die Funktionen mandelbrot und mandelbrot reell.)

3 Das Newtonfraktale – Pflichtaufgaben Teil 3

Das insbesondere in der numerischen Mathematik enorm wichtige Newton-Verfahren kann ebenfalls zur Erzeugung von Fraktalen benutzt werden. Bei dem Newton-Verfahren handelt es sich um einen grundlegenden Algorithmus zur Bestimmung von Lösungen der Gleichung f(x) = 0 für f : ℝ^{m} → ℝ^{m}. Das Verfahren besteht darin, dass man für einen geeigneten Startwert x_{0} ∈ ℝ^{m} die Iterationsvorschrift

x_{n+1} = x_{n} − Df(x_{n})^{−1}f(x_{n}), n = 0, 1, 2, ...

solange anwendet, bis der Defekt ∥f(x_{n})∥ unter einer vorgegebenen Toleranz liegt oder eine maximale Anzahl Iterationen erreicht ist. Der Index n bezeichnet hierbei den Iterationsindex und Df(x) die Jacobimatrix (\frac{∂f_{i}(x)}{∂x_{j}})_{i,j} an der Stelle x ∈ ℝ^{m}. Eine genauere Beschreibung findet sich in (vermutlich) jedem einführenden

Buch in die numerische Mathematik, siehe auch http://de.wikipedia.org/wiki/Newton Verfahren.

In Matlab kann man das mit Hilfe des Backslash Operators \ wie folgt implementieren:

[figure]

  tol     = 1e-8;          % Toleranz
  maxsteps = 8;            % max. Anzahl Iterationen
  x       = x0;            % Startwert der Iteration
  defekt  = norm(f(x));    % vorgegebener Defekt
  steps   = 0;             % Start im Nullten Schritt
  while defekt > tol && steps < maxsteps
    xx = x - Df(x) \ f(x);
    x = xx;                % x aufdatieren
    defekt = norm(f(x));   % Defekt berechnen
    steps = steps + 1;     % Schritte weiterzaehlen

  end

[figure]

Beachte: x = A\b berechnet die Lösung x der Matrixgleichung Ax = b.

Aufgabe 8. Erstelle aus obigem Code eine Funktion newton, die nach Eingabe der Funktionshandles für f, Df und des Startwertes x0 die Newtoniteration durchführt und das Ergebnis der Iteration x_{n} (also den letzten x-Wert), sowie die benötigte Anzahl Schritte zurückgibt. Wähle für die Parameter tol und maxiter des Verfahrens die obigen Werte.

Aufgabe 9. Probiere die Funktion newton an den folgenden Beispielen aus

f(x) = (x − 2)(x + 3), f(x) = sin(x)e^{x}.

Starte mit verschiedenen Werten x_{0} ∈ ℝ und beobachte das Verhalten des Algorithmus (zum Beispiel im Debugger). Zeichne auch den Graphen von f und lokalisiere die Nullstellen.

Aufgabe 10 (Pflichtaufgabe). Wende das Newtonverfahren auf die reelle Version der komplexen Funktion f : ℂ → ℂ, z ↦ z^{3} − 1 an. Diese ist gegeben durch

f(x, y) = (\begin{matrix}x^{3} − 3xy^{2} − 1 \\ 3x^{2}y − y^{3} \end{matrix}) .

Die Jacobimatrix Df(x, y) ist hier

Df(x, y) = (\begin{matrix}3(x^{2} − y^{2}) & −6xy \\ 6xy & 3(x^{2} − y^{2}) \end{matrix}) .

Starte hierbei mit verschiedenen Werten (x_{0}, y_{0}) in der Nähe der drei Nullstellen

(x_{1}, y_{1}) = (1, 0), (x_{2}, y_{2}) = (−\frac{1}{2}, \frac{\sqrt{3}}{2}), (x_{3}, y_{3}) = (−\frac{1}{2}, −\frac{\sqrt{3}}{2}).

Aufgabe 11 (Pflichtaufgabe). Auch zum Newtonverfahren gibt es ein Fraktal, das sogenannte Newton-Fraktal, das wie folgt erzeugt wird: Führe das obige Newtonverfahren für verschiedene Anfangswerte in einem Bereich, der die Nullstellen umfasst durch. Färbe die Anfangswerte (x_{0}, y_{0}) ∈ ℝ^{2} entweder danach,

• gegen welche der drei Nullstellen die Iteration konvergiert

oder

• nach wieviel Schritten die Iteration konvergiert (wie bei der Gefangenenmenge)

und plotte das Ergebnis (z.B. wie in Aufgabe 5 mit pcolor). Zeichne das Newton Fraktal für mindestens eins der folgenden Polynome

p_{1}(z) = (z + 1)(z − 1)(z + 1 + i)(z + 1 − i)  p_{2}(z) = z(z − 1)(z − i)(z − i + 1)
p_{3}(z) = (z + 1)(z − 1)(z + 1 + i)(z + 1 + i)(z + 1 − i)  p_{4}(z) = z^{5} + 4z.