Das 3k+1(Collatz)Problem

Wir betrachten folgendes Problem. Für jedes nϵℕ ist die Funktion f definiert durch

f:=(Piecewise[{{(n/2), falls n gerade ist}, {3n+1, sonst}}])

eine derartige Funktion ist in Mathematica einfach zu implementieren

collatz[n_]:=If[EvenQ[n],(n/2),3*n+1]

collatz[7]

Die Frage die wir nun stellen ist, was bei wiederholtem, iteriertem Anwenden der Funktion f passiert. Betrachtet man z.B. die Zahl 5 als Startzahl, so erhält man durch wiederholtes Anwenden der Funktion f die Folge der Zahlen 5,16,8,4,2,1 Also terminiert die Folge der Iterierten von f bei Eingabe von 5 mit dem Wert 1. Es ist bis heute nicht bewiesen, ob dies für jede beliebige Eingabe für die Funktion f auch gilt. Wir werden nun eine Mathematica Funktion schreiben, die für jeden Eingabewert die Folge der Iterierten bei Anwendung der Funktion f auf diese Eingabezahl ausgibt. Wir gehen dabei auf zwei unterschiedliche Arten vopr. Zunächst halten wir uns bei der Implementierung an einen Programmierstiel der sich am "klassichen" Programmierstil orientiert. Danach wird eine Modellierung des Problems gezeigt, die die speziellen (nicht nur mathematischen) Fähigkeiten von Mathematica verwendet.

Der klassische Ansatz

Wir beginnen mit der Iteration der Funktion Test, dabei wird als Maximalzahl der möglichen Iterationen1000 angenommen. Dazu definieren wir eine Funktion iteriereTest die genau dies tut, dabei verwenden wir eine klassische, kopfgesteuerte while-Schleife:

iteriereCollatz[n_]:=Module[{maxiter=1000,iter=1, out={n}, temp=n},
While[
     collatz[temp]>1 && iter<maxiter,
    {iter=iter+1,temp=collatz[temp],AppendTo[out,temp]} 
]
;If[iter<maxiter,AppendTo[out,1],];out
]

iteriereCollatz[7]

Length @ %

Die Anzahl der Iterationen, bis die Folge bei 1 anlangt erhalten wir nun einfach durch die Länge der Ausgabeliste der Funktion iteriereTest. Damit ist es einfach für einen gegebenen Abschnitt der natürlichen Zahlen eine Funktion zu definieren, die für jede Zahl, die Länge der iterierten Folge liefert

folgenLänge[min_,max_]:= Module[{out={}},
For[n=min,n<max,n++,AppendTo[out,{n,Length[iteriereCollatz[n]]}]];out  
]

folgenLänge[3,11]

Darauf aufbauen ist es nun sehr einfach eine Funktion zu definieren, die für ein gegebenes Intervall natürlicher Zahlen die Länge der iterierten Folgen des 3k+1-Problems graphisch darstellt, dazu dient die Funktion ListPlot

zeigeIterierte[min_,max_]:=Module[{n,liste={} },
For[n=min,n<max,n++,
AppendTo[ liste,{n,Length[iteriereCollatz[n]]} ]];ListPlot[liste, PlotStyle ->PointSize[0.005]]]

zeigeIterierte[3,1000]

Interessant ist dabei die Bildung "horizontaler Linien", die anzeigt, daß hier Zahlen nebeneinander, oder zumindest in der Nähe liegen, die die gleiche Anzahl von Iterationen benötigen, um die Folge zum Terminieren zu bringen. Deutlicher zu sehen ist die bei der Analyse größerer Zahlenintervalle

{zeit,plot} = zeigeIterierte[1000000,1010000]//Timing;

Sow[plot]

zeit

Remove[plot,zeit,zeigeIterierte,folgenLänge,collatz,iteriereCollatz]

Der Mathematica Ansatz

Die oben definierte Funktion iteriereTest löst das Problem der Iteration auf die klassische Weise im Sinne der prozeduralen Programmierung. Wir gehen nun den "Mathematica-Weg". Das Iterieren einer Funktion unter Beibehaltung aller Zwischenwerte wird mit der Funktion NestWhileList realisiert. Wir definieren zunächst die Funktion Collatz, die der Funktion testn enstpricht (Das 3k+1 Problem heißt, zu Ehren des Mathematikers Collatz, der es als erster formulierte, Collatz-Problem).

Hier der direkt Ansatz in Mathematica-Manier:

Der direkte Ansatz

collatz[n_]:=If[EvenQ[n],(n/2),3*n+1]

Die Funktion NestWhileList iteriert die Funktion Collatz so lange das Ergebnis der Funktoin Collatz größer als Eins ist:

NestWhileList[collatz,27,#>1&]

Nun muß nur noch die Funktion definiert werden, die für eine gegebene Zahl diese Collatz-Iteration durchführt:

collatzFolge[n_]:= NestWhileList[collatz,n,#1>1&]

collatzFolge[11]

Nun ist die Lösung des Problems (oder besser die graphische Darstellung desselben, denn gelöst ist es wie gesagt, noch nicht) mit Mathematica Bordmitteln einfach. Die Funktion CollatzFolge wird einfach auf die zu untersuchende Folge von Zahlen gemapt.

collatzFolge /@ Range[2,10]// TableForm

Die graphische Darstellung ist ebenso einfach. Vom letzen Ergebnis wird einfach jeweils die Folgenlänge ermittelt. Damit kann mit einfachen Mitteln sofort die jeweilige Folgenlänge von Collatzfolgen für alle Zahlen in einem gegebenen Intervall berechnet werden. Dies liefert schließlich die Eingabe für die Funktion ListPlot.

Length /@ (collatzFolge /@ Range[2,10])

Dieses Ergebnis liefert schließlich die Eingabe für die Funktion ListPlot, die das Ergebnis graphisch darstellt.  Die Mathematica-nahe Implementierung ist nicht nur kürzer, sondern, was die Laufzeit angeht, auch die bei weitem bessere Lösung

collatzPlot[n_,m_]:=ListPlot[ Length /@ (collatzFolge /@ Range[n,m]),PlotStyle ->PointSize[0.005]]/;n<m

{zeit,plot}= collatzPlot[1000000,1010000]//Timing;

Show[plot]

Neben dieser Beobachtung ist es auch interessant zu untersuchen, wie häufig bestimmte Folgenlängen im Rahmen einer Collatziteration eines bestimmten Intervalls auftreten. Hierzug definieren wir zunächst eine Funktion, die die Länge der iterierten Folge direkt ausgibt:

folgenLänge[n_]:= Length @ collatzFolge[n]

Damit ist es einfach, für ein gegebenes Intervall die jeweilige Länge der Collatzfolge zu bestimmen. Wir betrachten in diesem Falle Paare {n,l} wobei  n eine natürliche Zahl und l die zugehörige Länge der Colatzfolge ist:

folgenLänge[27]

collatzFolge[27]

temp=Table[{n,folgenLänge[n]},{n,100,110}]

Damit ist es einfach zu bestimmen, wie die Länge der Colatzfolgen im gegebenen Intervall verteilt ist (die Anweisung {x_,0}->Sequence[] bewirkt dabei, daß alle Einträge, die nicht als Länge einer Colatzfolge auftreten unterdrückt werden.):

sequenzen[n_,m_]:=Module[{min, max,liste,temp},
liste = Table[{i,folgenLänge[i]},{i,n,m}];
temp=Last[#]& /@ liste ;
{min,max}= {Min[temp],Max[temp]};
max=Last[#]& /@ liste // Max;
{Table[{i,Length @ Select[liste,Last[#]==i&]},{i,min,max}]/. {x_,0}->Sequence[],min,max}
]

temp =sequenzen[2,100000];

ListPlot[temp[[1]], PlotRange->All, AxesOrigin->{0,0},AxesLabel->{Style["Periodenlänge",Larger, Bold], Style["Häufigkeit",Larger,Bold]}]

Auch hier zeigt sich eine sehr interessante Verteilung der Häufigkeit der einzelnen Folgenlängen in der Ausprägung zweier Maxima, eines zwischen 50 und 80 ein zweites zwischen 140 und 160. Diese Ausprägung (zweier) unterschiedlicher Maxima ist auch auf einzelnen Teilintervallen zu beobachten:

plotliste=Table[sequenzen[i,i+1000][[1]],{i,2,10000,1000}];

ListPlot /@ plotliste

resp:

ListPlot[plotliste, Background->LightGray]

Auch die Entwicklung der Folgenlänge für ausgewählte Sequenzen ist einer Betrachtung wert. Im Folgenden ist die Entwicklung der Folgenlänge für (3^n), (5^n)  und (7^n) dargestellt. Auffällig ist hier die “Bandbildung” bei der Folgenlänge der Folge (3^n)

temp=Map[folgenLänge,Table[{3,5,7}^n,{n,1,500}]\[Transpose],{2}];

ListPlot[temp, Background->Gray,PlotStyle->{Red, Blue,Cyan}]

Der Nachteil bei dieser Lösung ist, daß in der Rekursion viele Berechnungen doppelt durchgeführt werden. Wenn etwa collatz[16] aufgerufen wird, so werden die erforderlcihen Berechnungen alle ausgeführt, obwohl diese möglicherweise schon bei einer früherern Bwerechnung aufgetreten sind. Das läßt sich vermeiden

Vermeiden doppelter Berechnungen

Statt wie oben definieren wir:

collatzRek[n_]:=collatzRek[n]=If[EvenQ[n],n/2,3*n+1]

und

collatzSequence[n_]:=NestWhileList[collatzRek,n,#>1&]

sowie

collatzLength[n_]:=Length@collatzSequence[n]

Damit geht es etwas flotter:

seq=Table[2^i+1,{i,1,1000}];First @Timing@(folgenLänge[#]&/@ seq)

First @ Timing@(collatzLength[#]&/@ seq)