Zufallswanderer Modelle

Simulation komplexer Systeme

H.J. Lüdde

CA Bibliothek

> restart:

> with(plots):

> with(plottools):

> with(linalg):

Zweifarbige Darstellung des CA Feldes

> plotCA2 := proc(A,m,n)

> global sumCA;

> local i,j,ij,k,r,CA,board;

> CA:=array(1..m*n):

> r := plottools[rectangle]([0,1],[1,0],color=blue):

> # Zeichne weiße Grundfläche

> board:=plottools[rectangle]([0,m],[n,0],color=white):

> # Zeichne farbige Zellen

> ij:=0:

> for i from 0 to m-1 do

> for j from 0 to n-1 do

> if A[m-i,j+1]>0 then

> ij:=ij+1:

> CA[ij]:=plottools[translate](r,j,i):

> fi:

> od:

> od:

> sumCA:=ij:

> # Zeichne CA Feld

> if ij=0 then

> board:

> else

> seq(CA[k],k=1..ij),board:

> fi;

> end:

Dreifarbige Darstellung des CA Feldes

> plotCA3 := proc(A,m,n)

> global sumCA;

> local i,j,ij,k,r,CA,board;

> CA:=array(1..m*n):r:=array(1..2):

> r[1] := plottools[rectangle]([0,1],[1,0],color=green):

> r[2] := plottools[rectangle]([0,1],[1,0],color=red):

> # Zeichne weiße Grundfläche

> board:=plottools[rectangle]([0,m],[n,0],color=white):

> # Zeichne farbige Zellen

> ij:=0:

> for i from 0 to m-1 do

> for j from 0 to n-1 do

> if A[m-i,j+1]=1 then

> ij:=ij+1:

> CA[ij]:=plottools[translate](r[1],j,i):

> elif A[m-i,j+1]>1 then

> ij:=ij+1:

> CA[ij]:=plottools[translate](r[2],j,i):

> fi:

> od:

> od:

> sumCA:=ij:

> # Zeichne CA Feld

> if ij=0 then

> board:

> else

> seq(CA[k],k=1..ij),board:

> fi;

> end:

Prozedur für binäre Zufallszahlen

> oz := proc() rand() mod 2 end:

Prozedur für binäre Zufallszahlen mit unterschiedlicher Dichte

(Verteilung: (100/k)% Felder werden mit 1 belegt)

> oz1 := proc(k) iquo((rand() mod k)+1,k) end:

Prozedur für Zufallszahlen aus [0,1]

> oz01 := proc() rand()/10^12 end:

Besetzung einer mxn Matrix mit Zufallszahlen

(k ist die Dichte der Verteilung: (100/k)%)

> RandCA:= proc(m,n,k)

> local i,j;

> for i from 1 to m do

> row(i):=seq(oz1(k),j=1..n):

> od:

> matrix([seq([row(j)],j=1..m)]);

> end:

Regeln für mehrfache Zufallswanderer

> RuleMWalker:= proc(A,pA,i,j,m,n)

> local iu,iuu,ju,juu,io,ioo,jo,joo,a,die;

> # Periodische Randbedingungen

> if i=1 then iu:=m:iuu:=m-1:

> elif i=2 then iu:=1:iuu:=m:

> else iu:=i-1:iuu:=i-2:

> fi:

> if j=1 then ju:=n:juu:=n-1:

> elif j=2 then ju:=1:juu:=n:

> else ju:=j-1:juu:=j-2:

> fi:

> if i=m then io:=1:ioo:=2:

> elif i=m-1 then io:=m:ioo:=1:

> else io:=i+1:ioo:=i+2:

> fi:

> if j=n then jo:=1:joo:=2:

> elif j=n-1 then jo:=n:joo:=1:

> else jo:=j+1:joo:=j+2:

> fi:

> # 'Walker' Regel in 'v. Neumann'- Umgebung

> die:=rand(2..5):

> a:=A[i,j]:

> if A[i,j]=0 then

> if A[io,j]=2 or A[i,ju]=3 or A[iu,j]=4 or A[i,jo]=5 then

> a:=die():

> else

> a:=0:

> fi:

> elif A[i,j]=2 then # Wanderer will in Zelle A[iu,j]

> if A[iu,j]>0 then # Zielzelle besetzt

> a:=die():

> else

> if (A[iu,ju]=3 and pA[iu,ju]>pA[i,j]) or (A[iuu,j]=4 and pA[iuu,j]>pA[i,j]) or

> (A[iu,jo]=5 and pA[iu,jo]>pA[i,j]) then

> a:=die():

> else

> a:=0:

> fi:

> fi:

> elif A[i,j]=3 then # Wanderer will in Zelle A[i,jo]

> if A[i,jo]>0 then # Zielzelle besetzt

> a:=die():

> else

> if (A[iu,jo]=4 and pA[iu,jo]>pA[i,j]) or (A[i,joo]=5 and pA[i,joo]>pA[i,j]) or

> (A[io,jo]=2 and pA[io,jo]>pA[i,j]) then

> a:=die():

> else

> a:=0:

> fi:

> fi:

> elif A[i,j]=4 then # Wanderer will in Zelle A[io,j]

> if A[io,j]>0 then # Zielzelle besetzt

> a:=die():

> else

> if (A[io,jo]=5 and pA[io,jo]>pA[i,j]) or (A[ioo,j]=2 and pA[ioo,j]>pA[i,j]) or

> (A[io,ju]=3 and pA[io,ju]>pA[i,j]) then

> a:=die():

> else

> a:=0:

> fi:

> fi:

> elif A[i,j]=5 then # Wanderer will in Zelle A[i,ju]

> if A[i,ju]>0 then # Zielzelle besetzt

> a:=die():

> else

> if (A[io,ju]=2 and pA[io,ju]>pA[i,j]) or (A[i,juu]=3 and pA[i,juu]>pA[i,j]) or

> (A[iu,ju]=4 and pA[iu,ju]>pA[i,j]) then

> a:=die():

> else

> a:=0:

> fi:

> fi:

> fi:

> a:

> end:

Prioritätenbesetzung

> priority:=proc(A,i,j,m,n)

> if A[i,j]>0 then oz01(): else 0: fi:

> end:

Warning, new definition for norm

Warning, new definition for trace

>

1) Der einzelne Zufallswanderer

Der Zufallswanderer startet in der Mitte des CA Feldes und kann sich in die vier Richtungen der Neumann Umgebung bewegen. Die Richtung wird mit Hilfe der Prozedur 'die' gewürfelt. Die erlaubten Zellenzustände sind: 0 (weiß): die Zelle ist leer und war noch nie besetzt, 1 (grün): die Zelle ist zur Zeit leer, war aber in der Vergangenheit schon mindestens einmal besetzt, 2-5 (rot): die Zelle ist besetzt und der Wanderer wird im nächsten Schritt in die Richtungen N(2), O(3), S(4), W(5) gehen. An den Rändern des CA Feldes gelten periodische Randbedingungen.

Regeln für einzelnen Zufallswanderer

> RuleSWalker:= proc(A,i,j,m,n)

> local iu,ju,io,jo,a,die;

> # Periodische Randbedingungen

> if i-1=0 then iu:=m; else iu:=i-1; fi;

> if j-1=0 then ju:=n; else ju:=j-1; fi;

> if i+1>m then io:=1; else io:=i+1; fi;

> if j+1>n then jo:=1; else jo:=j+1; fi;

> # 'Walker' Regel in 'v. Neumann'- Umgebung

> die:=rand(2..5):

> a:=A[i,j]:

> if A[i,j]<2 then

> # if A[io,j]=2 then a:=die(): fi:

> if A[io,j]=2 or A[i,ju]=3 or A[iu,j]=4 or A[i,jo]=5 then a:=die(): fi:

> else

> a:=1:

> fi:

> a:

> end:

Start des Zufallswanderers in der Mitte des CA Feldes

> InitSWalker:= proc(m,n)

> local B,die;

> die:=rand(2..5):

> B:=matrix(m,n,0): # Initialisierung der unbesetzten Zellen

> B[iquo(m,2)+1,iquo(n,2)+1]:=die(): # Wahl einer Zufallsrichtung des Wanderers

> B:

> end:

Beginn des Programmes

1) Initialisierung

> M:=21: # Zeilendimension

> N:=21: # Spaltendimension

> A:=InitSWalker(M,N): # Vorgabe des CA für Zufallswanderer

2) Propagation über nstep Zeitschritte

> nstep:=100:

> p.1:=plotCA3(A,M,N):

> for k from 2 to nstep do

> for i from 1 to M do

> row(i):=seq(RuleSWalker(A,i,j,M,N),j=1..N):

> od:

> A:=matrix(M,N,[seq(row(j),j=1..M)]):

> p.k:=plotCA3(A,M,N):

> od:

> k:=k-1:

> display(seq(display(p.j),j=1..k),insequence=true,scaling=constrained,axes=none); # Animation der Propagation

Aktiviere die Grafik durch 'anklicken' und starte die Animation. Die rote Zelle zeigt die aktuelle Position, die grünen Flächen die Spur des in der Mitte des CA gestarteten Zufallswanderers.

2) Diffusionsautomat

Mit Hilfe dieses Automaten wird die Diffusion einer zunächst in der Mitte des CA Feldes konzentrierten Substanz modelliert. Das erreicht man durch die gleichzeitige Propagation von mi x ni Zufallswanderern. Dabei ist zu beachten, dass nicht zur selben Zeit zwei Wanderer die selbe Zelle besetzen (die Anzahl der Wanderer muss erhalten bleiben) dürfen. Deshalb wird jedem Wanderer eine Priorität zugeordnet (Zufallszahl [0..1]). Wollen zwei Wanderer die selbe leere Zelle besetzen, so hat derjenige mit der höheren Priorität den Vorrang.

Start des Diffusionsautomaten in der Mitte des CA Feldes

> InitMWalker:= proc(m,n,mi,ni)

> local B,die,dm,dn,i,j;

> die:=rand(2..5):

> dm:=iquo(m-mi,2):

> dn:=iquo(n-ni,2):

> B:=matrix(m,n,0): # Initialisierung der unbesetzten Zellen

> for i from dm+1 to dm+mi do

> for j from dn+1 to dn+ni do

> B[i,j]:=die(): # Wahl einer Zufallsrichtung des Wanderers

> od:

> od:

> B:

> end:

Initialisierung der Prioritätenmatrix

> InitPriority:= proc(m,n,mi,ni)

> local pB,dm,dn,i,j;

> dm:=iquo(m-mi,2):

> dn:=iquo(n-ni,2):

> pB:=matrix(m,n,0): # Initialisierung der unbesetzten Zellen

> for i from dm+1 to dm+mi do

> for j from dn+1 to dn+ni do

> pB[i,j]:=oz01(): # Zufalsswert aus [0,1]

> od:

> od:

> pB:

> end:

Beginn des Programmes

1) Initialisierung

> M:=40:MI:=10: # Zeilendimension

> N:=40:NI:=10: # Spaltendimension

> A:=InitMWalker(M,N,MI,NI): # Vorgabe des CA für Diffusionsautomat

> pA:=InitPriority(M,N,MI,NI): # Vorgabe der Prioritäten der besetzten Zellen

> display(plotCA2(A,M,N),scaling=constrained,axes=none);# Grafische Darstellung der Anfangsbedingung

2) Propagation über nstep Zeitschritte

> nstep:=200:

> p.1:=plotCA2(A,M,N):

> for k from 2 to nstep while sumCA>0 do

> for i from 1 to M do

> rowA(i):=seq(RuleMWalker(A,pA,i,j,M,N),j=1..N):

> od:

> A:=matrix(M,N,[seq(rowA(j),j=1..M)]):

> for i from 1 to M do

> rowpA(i):=seq(priority(A,i,j,M,N),j=1..N):

> od:

> pA:=matrix(M,N,[seq(rowpA(j),j=1..M)]):

> p.k:=plotCA2(A,M,N):

> od:

> k:=k-1:

> display(seq(display(p.j),j=1..k),insequence=true,scaling=constrained,axes=none); # Animation der Propagation

>

3) Aggregationsautomat (Diffusionsbegrenztes Wachstum)

Unter dem Begriff diffusionsbegrenztes Wachstum (Diffusion Limited Aggregation (DLA)) versteht man die Entstehung fraktaler Aggregate. In einem Medium lagern sich frei bewegliche Teilchen an einen wachsenden Cluster an. Dabei wird die Struktur des Wachstums vom Zufall bestimmt: die durch die Brown'sche Molekularbewegung wandernden Teilchen bleiben am ersten zufällig erreichten Kontaktpunkt mit dem Wachstumskeim hängen. Es bilden sich Aggregate fraktaler Dimension, einem wichtigen Phänomen natürlicher Strukturbildung. Das DLA Modell wurde 1981 von Witten und Sander (Phys. Rev. Lett. 47 , 1400) eingeführt.

Falls Bindung in Moore Umgebung zu Konsationskeim, setze A[i,j]=1

> aggregation:= proc(A,i,j,m,n)

> local a,iu,io,ju,jo:

> if i-1=0 then iu:=m; else iu:=i-1; fi;

> if i+1>m then io:=1; else io:=i+1; fi;

> if j-1=0 then ju:=n; else ju:=j-1; fi;

> if j+1>n then jo:=1; else jo:=j+1; fi;

> a:=A[i,j]:

> if A[i,j]>1 then

> if A[iu,ju]=1 or A[iu,j]=1 or A[iu,jo]=1 or

> A[i,ju]=1 or A[i,j]=1 or A[i,jo]=1 or

> A[io,ju]=1 or A[io,j]=1 or A[io,jo]=1 then a:=1: fi:

> fi:

> a:

> end:

Würfele neue Richtungen für Wanderer

> NewOrientation:= proc(A,i,j)

> local die,a:

> die:=rand(2..5):

> if A[i,j]>0 then a:=die(): else a:=A[i,j]: fi:

> a:

> end:

Beginn des Programmes

1) Initialisierung

> M:=51: # Zeilendimension

> N:=51: # Spaltendimension

> B:=RandCA(M,N,6): # Zufallsbelegung des CA mit 1 oder 0

> for i from 1 to M do # Falls A[i,j]=1 würfele Richtung

> rowA(i):=seq(NewOrientation(B,i,j),j=1..N):

> od:

> A:=matrix(M,N,[seq(rowA(j),j=1..M)]):

> for i from 1 to M do

> rowpA(i):=seq(priority(A,i,j,M,N),j=1..N):

> od:

> pA:=matrix(M,N,[seq(rowpA(j),j=1..M)]): # Setze Prioritätenmatrix

> A[iquo(M,2)+1,iquo(N,2)+1]:=1: # Kondensationskeim in der Mitte des CA

> display(plotCA3(A,M,N),scaling=constrained,axes=none);# Grafische Darstellung der Anfangsbedingung

2) Propagation über nstep Zeitschritte

> nstep:=100:

> p.1:=plotCA3(A,M,N):

> for k from 2 to nstep do

> for i from 1 to M do

> rowA(i):=seq(RuleMWalker(A,pA,i,j,M,N),j=1..N): # Berechne neue Positionen

> od:

> A:=matrix(M,N,[seq(rowA(j),j=1..M)]):#print(A);

> for i from 1 to M do

> rowB(i):=seq(aggregation(A,i,j,M,N),j=1..N): # Überprüfe Aggregationsbedingung

> od:

> A:=matrix(M,N,[seq(rowB(j),j=1..M)]):#print(A):

> for i from 1 to M do

> rowpA(i):=seq(priority(A,i,j,M,N),j=1..N): # Berechne Prioritätsmatrix

> od:

> pA:=matrix(M,N,[seq(rowpA(j),j=1..M)]):#print(pA);

> p.k:=plotCA3(A,M,N):#print(A);

> od:

> k:=k-1:

> display(seq(display(p.j),j=1..k),insequence=true,scaling=constrained,axes=none); # Animation der Propagation

>

>

4) Kondensationsautomat

Unter Kondensation versteht man den Übergang eines Stoffes in einen Aggregatzustand niedrigerer Gesamtenergie. Z.B. führt die Kondensation in einem mit Wasserdampf übersättigten Volumen zur Tröpfchenbildung. Der Kondensationsprozess hält solange an, bis die Gesamtenergie des Systems bestehend aus isolierten Molekülen (Monomeren) und Clustern (Anlagerung von zwei (Dimer) oder mehr Molekülen) ein Minimum erreicht. In der nachfolgenden Prozedur 'energy' wird die Gesamtenergie berechnet aus der Anzahl der vorhandenen Bindungen (Bindungsenergie ist proportional der Zahl der besetzten Nachbarzellen in Neumann Geometrie) abzüglich der Oberflächenenergie (Anzahl der unbesetzten Nachbarzellen). In jedem Zeitschritt wird die Wanderungsrichtung der Cluster gewürfelt, die Cluster entsprechend verschoben und die Gesamtenergie des Systems berechnet. Ist die Gesamtenergie kleiner als die bisher gespeicherte kleinste Energie, wird der neue Zustand übernommen, ansonsten verworfen.

Energie des Gesamtsystems (Tröpfchenmodell)

> energy:= proc(A,m,n,alpha,beta)

> local b,s,i,iu,io,j,ju,jo:

> b:=0:s:=0:

> for i from 1 to m do

> if i-1=0 then iu:=m; else iu:=i-1; fi;

> if i+1>m then io:=1; else io:=i+1; fi;

> for j from 1 to n do

> if j-1=0 then ju:=n; else ju:=j-1; fi;

> if j+1>n then jo:=1; else jo:=j+1; fi;

> if A[i,j]>0 then

> if A[iu,j]>0 then b:=b+1: else s:=s+1: fi:

> if A[i,jo]>0 then b:=b+1: else s:=s+1: fi:

> if A[io,j]>0 then b:=b+1: else s:=s+1: fi:

> if A[i,ju]>0 then b:=b+1: else s:=s+1: fi:

> fi:

> od:

> od:

> (-beta*b+(1-beta)*s)/alpha:

> end:

Würfele neue Richtungen für Wanderer

> NewOrientation:= proc(A,i,j)

> local die,a:

> die:=rand(2..5):

> if A[i,j]>0 then a:=die(): else a:=0: fi:

> a:

> end:

Beginn des Programmes

1) Initialisierung

> M:=20: # Zeilendimension

> N:=20: # Spaltendimension

> B:=RandCA(M,N,20): # Zufallsbelegung des CA mit 1 oder 0

> for i from 1 to M do # Falls A[i,j]=1 würfele Richtung

> rowA(i):=seq(NewOrientation(B,i,j),j=1..N):

> od:

> A:=matrix(M,N,[seq(rowA(j),j=1..M)]):

> for i from 1 to M do

> rowpA(i):=seq(priority(A,i,j,M,N),j=1..N):

> od:

> pA:=matrix(M,N,[seq(rowpA(j),j=1..M)]):

> display(plotCA2(A,M,N),scaling=constrained,axes=none);# Grafische Darstellung der Anfangsbedingung

2) Propagation über nstep Zeitschritte

> nstep:=50:

> alpha:=0.15:beta:=0.5:

> eold:=energy(A,M,N,alpha,beta):print(E=eold);

> p.1:=plotCA2(A,M,N):

> for k from 2 to nstep do

> for l from 1 to 1000 do

> for i from 1 to M do

> rowB(i):=seq(NewOrientation(A,i,j),j=1..N): # Würfele neue Orientierung der alten Positionen

> rowA(i):=seq(RuleMWalker(A,pA,i,j,M,N),j=1..N): # Berechne neue Positionen

> od:

> A:=matrix(M,N,[seq(rowA(j),j=1..M)]):

> enew:=energy(A,M,N,alpha,beta):#print(enew); # Berechne Energie bezogen auf neue Positionen

> if eold<= enew then A:=matrix(M,N,[seq(rowB(j),j=1..M)]): fi:#print(A); # Setze Matrix zurück

> for i from 1 to M do

> rowpA(i):=seq(priority(A,i,j,M,N),j=1..N):

> od:

> pA:=matrix(M,N,[seq(rowpA(j),j=1..M)]):#print(pA);

> if eold>enew then eold:=enew: break fi:

> od:

> p.k:=plotCA2(A,M,N):print(l, E=eold);

> if l=1001 then break fi:

> od:

> k:=k-1:

> display(seq(display(p.j),j=1..k),insequence=true,scaling=constrained,axes=none); # Animation der Propagation

>