SIMULATIONEN



 

Monte Carlo Simulation
 

 

Das klassische Beispiel zur Monte Carlo - Methode berechnet die Zahl Pi aus der Zahl der Zufallstropfen, die auf einen Viertelkreis fallen zur Gesamtzahl Tropfen, die man auf eine quadratische Fläche fallen lässt. Mit der Monte Carlo-Methode lassen sich aber auch beliebige Flächen berechnen.
 

Auf eine frei gewählte Polygonfläche lässt man eine grosse Zahl von Zufallspunkte fallen und zählt, wie viele davon auf die Polygonfläche fallen. Aus dem Verhältnis der Zahl der Punkte im Polygon zur Anzahl aller Punkte lässt sie die Polygonfläche berechnen.

Die Eckpunkte des Polygons werden mit dem linken Mausklick gewählt. Mit dem rechten Mausklick wird die Polygonfläche gefärbt und der Zufallsregen gestartet. Auf Grund der Pixelfarbe des Hintergrunds wird entschieden, ob der Punkt auf die Polygonfläche fällt.

  

Programm: [►]

# SimEx1.py

from gpanel import *
import random
import thread

NB_DROPS = 10000

def onMousePressed(x, y):
    if isLeftMouseButton():
        pt = [x, y]
        move(pt)
        circle(0.5)
        corners.append(pt)
    if isRightMouseButton():
        thread.start_new_thread(go,())

def go():
    global nbHit
    setColor("gray")
    fillPolygon(corners)
    title("Working. Please wait...")
    for i in range(NB_DROPS):
        x = 100 * random.random()
        y = 100 * random.random()
        color = getPixelColorStr(x, y)
        if color == "black":
            setColor("green")
            point(x, y)
        if color == "gray" or color == "red":
            nbHit += 1
            setColor("red")
            point(x, y)
    title("All done. #hits: " + str(nbHit) + " of " + str(NB_DROPS))

makeGPanel(0, 100, 0, 100, mousePressed = onMousePressed)
title("Select corners with left button. Start dropping with right button")
bgColor("black")
setColor("gray")
corners = []
nbHit = 0
keep()
Progammcode markieren (Ctrl+C copy, Ctrl+V einfügen)

  ► Erklärungen zum Programmcode
 

Mit isLeftMouseButton() bzw. isRightMouseButton() kann man herausfinden, ob der Mausevent mit der linken oder rechte Maustaste ausgelöst wurde.

 

Simulation von physikalischen Experimenten

 

 

In der Physik sind nummerische Simulationen oft einfacher als die mathematischen Herleitungen. Ausgehend vom aktuellen Zustand wird in regelmässigen Zeitabständen jeweils der neue Zustand berechnet.
 

In diesem Beispiel wird der vertikale Wurf simuliert. Ausgehend von der Startgeschwindigkeit v und der Beschleunigung a = -g mit g = 9.81, wird in Zeitintervallen dt aus der aktuellen Geschwindigkeit die neue Geschwindigkeit
v = v - g * dt
und aus dem aktuellen Ort der neue Ort
y = y + v * dt
berechnet.

 

  

Programm: [►]

# SimEx2.py

from gpanel import *
import time


g = 9.81
dt = 0.1 

makeGPanel(-1, 11, -220, 70)
drawGrid(0, 10, -200, 50, "gray")
setColor("red")
lineWidth(2)

t = 0; y = 0; v = 20
pos(t, y)
while t < 10:
	draw(t, y)
	v = v - g * dt
	y = y + v * dt	
	t = t + dt
	time.sleep(dt)
keep()
Progammcode markieren (Ctrl+C copy, Ctrl+V einfügen)

 

  ► Erklärungen zum Programmcode
 

time.sleep(dt) : Hält die Programmausführung während 0.1 s an. Es ergibt sich damit nur ungefähr eine Periode von 100 ms, da die Zeit für das Durchlaufen der anderen Anweisungen des Wiederholblocks nicht berücksichtigt wird.

Mit pos(t, y) wird der Grafikcursor an den ersten Punkt der Grafik gesetzt. Die nachfolgenden Werte werden mit draw(t, y) verbunden.

 

 

Gedämpfte harmonische Schwingung
 

Mit dem Newtonschen Gesetz F= ma berechhnet man aus der bekannten Kraft die Beschleunigung und verwendet wie vorher die Linearsierung.

Eine wesentliche Verbesserung der Konvergenz erreicht man, indem man für das erste Zeitintervall nur den halben Zeitschritt wählt (Halbschrittbeginn)

dt = 0.01
m = 0.5
k = 4
r = 0.1
t = 0; y = 0.8; v = 0
F = -k·y - r·v
a = F/m
v = v + a·dt
y = y + v·dt
 Zeitintervall in s
Masse
Federkonstante
Reibungskoeffizient
Anfangsbedingungen
Kraft
Beschleunigung
Neue Geschwindigkeit
Neue Koordinate
 

Programm: [►]

# SimEx3.py

from gpanel import *
import time

dt = 0.01
m = 0.5
k = 4
r = 0.1

t = 0; y = 0.8; v = 0

makeGPanel(-1, 11, -1.2, 1.2)
drawGrid(0, 10.0, -1.0, 1, "darkgray")
setColor("red")
lineWidth(3)

pos(t, y)
while t < 10:
      pos(t, y) if t == 0 else draw(t, y)
      F = -k*y - r*v
      a = F/m
      if t == 0:
          v = v + a*dt/2
      else:
          v = v + a*dt
      y = y + v*dt
      t = t + dt
      time.sleep(dt)
keep()
Progammcode markieren (Ctrl+C copy, Ctrl+V einfügen)

 

 

Feigenbaum-Diagramm
 

Für verschiedene Werte von r > 0 wird die die Konvergenz der Folge

an+1 = r * an * (1 - an) mit a0 = 0.5

untersucht.

In Programm wird jeweils aus der aktuellen Grösse x der neue Wert berechnet:

xnew = r * x * (1-x)

Für r < 1 gibt es einen Häufungspunkt  bei 0, die Folge konvergiert also gehen 0. Im Bereich zwischen 1 und 3  konvergiert die Folge ebenfalls. Für noch grössere r gibt es vorerst zwei und später mehr Häufungspunkte, die Folge konvergiert also nicht mehr. Für noch grössere Werte von r springt die Folge chaotisch hin und her.

 

Programm: [►]

# Feigenbaum.py

from gpanel import *

def f(x, r):
     return r * x * (1 - x)

makeGPanel(-0.6, 4.4, -0.1, 1.1)
title("Tree Diagram")
drawGrid(0, 4.0, 0, 1.0, "gray")
enableRepaint(False)
for z in range(1001):
    r = 4 * z / 1000
    a = 0.5
    for i in range(1001):
        a = f(a, r)
        if i > 500:
            point(r, a)
    repaint()
keep()
Progammcode markieren (Ctrl+C copy, Ctrl+V einfügen)

 

  ► Erklärungen zum Programmcode
 

enableRepaint(False): Deaktiviert das automatische Rendernbei jeder Grafikoperation. Mit repaint() führt man das Rendern dann an geeignter Stelle "von Hand" durch.. Dadurch wird die Grafik wesentlich schneller dargestellt da nicht bei jeder Grafikoperation (hier dem Zeichnen eines einzelnen Punkts) der ganze Bildschirm neu aufgebaut wird.