{ "cells": [ { "cell_type": "markdown", "id": "c9115fd9", "metadata": {}, "source": [ "\n", "\n", "\n", "\n", "\n", "\n", "
Institut für Theoretische Physik
Universität zu Köln 		 Prof. Dr. Simon Trebst
Christoph Berke
\n", "
\n", "

Statistische Physik

\n", "

Übungsblatt 3

\n", "
\n", "

Wintersemester 23/24

\n", "\n", "\n", "**Website:** https://www.thp.uni-koeln.de/trebst/Lectures/2023-StatPhys.shtml\n", "\n", "**Abgabe**: 30.10.2023, 10:00 Uhr **Besprechung**: 31.10.2023 \n", "\n", "**Name**: Bitte geben Sie Ihren Namen an. " ] }, { "cell_type": "markdown", "id": "73c78fa0", "metadata": {}, "source": [ "\n", "\n", "
Aufgabe 13: Simulation eines zweidimensionalen Gases, Teil 2 [4 Punkte]
" ] }, { "cell_type": "markdown", "id": "778d53b3", "metadata": {}, "source": [ "Wie in der vergangenen Woche finden Sie nachfolgend alle Funktionen implementiert, die für die Simulation eines zweidimensionalen Gases unter Berücksichtigung von elastischen Kollisionen zwischen den Teilchen benötigt werden. Führen Sie die nachfolgenden Zellen aus und bearbeiten Sie anschließend die unten stehende, ausführliche Aufgabenstellung." ] }, { "cell_type": "code", "execution_count": null, "id": "ba5a89c9", "metadata": {}, "outputs": [], "source": [ "using Random, LinearAlgebra, CairoMakie" ] }, { "cell_type": "code", "execution_count": null, "id": "a2c40f09", "metadata": {}, "outputs": [], "source": [ "# Initialize positions and velocities of particles in a 2D box.\n", "# x, y = Position arrays,\n", "# vx, vy = velocity arrays,\n", "# r = particle radius array\n", "# v = particle velocity (abs. value)\n", "# box = tuple with boundaries\n", "function initialize_particles!(x, y, vx, vy, r, v, box)\n", " \n", " # extract boundaries\n", " xmin, xmax, ymin, ymax = box\n", " \n", " # initialize each particle individually\n", " for i in 1:length(x)\n", " \n", " # random position with distance r to boundary\n", " x[i] = xmin+2*r[i] + rand() * (xmax - xmin - 4*r[i])\n", " y[i] = ymin+2*r[i] + rand() * (ymax - ymin - 4*r[i])\n", " \n", " # random velocity direction\n", " angle = rand()*2*pi\n", " vx[i], vy[i] = cos(angle)*v, sin(angle)*v\n", " \n", " end\n", " \n", "end\n", "\n", "# Same as above, but take account of the finite dimension of each particle. \n", "function initialize_particles_without_overlap!(x, y, vx, vy, r, v, box)\n", " \n", " # extract boundaries\n", " xmin, xmax, ymin, ymax = box\n", " \n", " # initialize each particle individually\n", " for i in 1:length(x)\n", " \n", " # reset counter of maximal tries\n", " tries, maxtries = 0, 1000\n", " \n", " # random position with distance r to boundary and no overlap to other particles\n", " overlapping = true\n", " while overlapping\n", " x[i] = xmin+2*r[i] + rand()*(xmax - xmin - 4*r[i])\n", " y[i] = ymin+2*r[i] + rand()*(ymax - ymin - 4*r[i])\n", " overlapping = false\n", " for j in 1:i-1\n", " if (x[i]-x[j])^2 + (y[i]-y[j])^2 < 1.01*(r[i]+r[j])^2\n", " overlapping = true\n", " end\n", " end\n", " tries += 1\n", " # maybe break loop because of maximum tries\n", " if tries > maxtries\n", " @error \"Could not place particle! Exceeded $(max_tries) tries\"\n", " return\n", " end\n", " end\n", " \n", " # random velocity direction\n", " angle = rand()*2*pi\n", " vx[i], vy[i] = cos(angle)*v, sin(angle)*v\n", " \n", " end\n", " \n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "b8a14600", "metadata": {}, "outputs": [], "source": [ "function move_particles!(x, y, vx, vy, dt, r, box)\n", " \n", " # move all particles by dt\n", " x .+= vx .* dt\n", " y .+= vy .* dt\n", " \n", " # check collision with bounding box\n", " for i in 1:length(x)\n", " # left / right\n", " if x[i]-r[i] < box[1]\n", " vx[i] *= -1\n", " elseif x[i]+r[i] > box[2]\n", " vx[i] *= -1\n", " end\n", " # up / down\n", " if y[i]-r[i] < box[3]\n", " vy[i] *= -1\n", " elseif y[i]+r[i] > box[4]\n", " vy[i] *= -1\n", " end\n", " end\n", " return 0\n", "end" ] }, { "cell_type": "code", "execution_count": null, "id": "2a856eb1", "metadata": {}, "outputs": [], "source": [ "function collide_particles!(x, y, vx, vy, m, r)\n", " \n", " # check all particle pairs (i,j) for collision\n", " for i in 1:length(x)\n", " for j in 1:i-1\n", " # particles i and j are too close --> collision\n", " if (x[i]-x[j])^2 + (y[i]-y[j])^2 < (r[i]+r[j])^2\n", " er = [x[j]-x[i], y[j]-y[i]]\n", " er = er ./ norm(er)\n", " vi = er[1]*vx[i] + er[2]*vy[i]\n", " vj = er[1]*vx[j] + er[2]*vy[j]\n", " vip = 2*(m[i]*vi + m[j]*vj)/(m[i]+m[j]) - vi\n", " vjp = 2*(m[i]*vi + m[j]*vj)/(m[i]+m[j]) - vj\n", " vx[i] = vx[i] + (vip - vi)*er[1]\n", " vy[i] = vy[i] + (vip - vi)*er[2]\n", " vx[j] = vx[j] + (vjp - vj)*er[1]\n", " vy[j] = vy[j] + (vjp - vj)*er[2]\n", " end\n", " \n", " end\n", " end\n", " \n", "end" ] }, { "cell_type": "markdown", "id": "de83d77f", "metadata": {}, "source": [ "### Vorbereitung: Aufenthaltswahrscheinlichkeit eines einzelnen Teilchens. ###\n", "Wir betrachten jetzt die Aufenthaltswahrscheinlichkeit eines **einzelnen** Teilchens in der Box. Dazu unterteilen wir (in Gedanken) das Gesamtvolument in einen linken und rechten Teil mit den jeweiligen Volumina $p V$ und $(1-p) V$ (Achtung: $p$ ist nicht der Druck. Später wird klar, warum dieser Parameter trotzdem $p$ heißt). In der folgenden Zelle betrachten wir mehrere, jeweils fixierte Werte von $p$. Wir lassen die Simulation jeweils für eine gewisse Zeitspanne ```twait``` laufen, und messen anschließend, ob sich das Teilchen im linken oder rechten Teilvolumen befindet. Diese Messung wiederholen wir mehrfach (```Nmeasure``` mal). Anschließend stellen wir die Häufigkeiten, mit denen das Teilchen im linken Teilvolumen gemessen wird als Funktion der Volumengröße $p$ dar." ] }, { "cell_type": "code", "execution_count": null, "id": "30cce5bf", "metadata": {}, "outputs": [], "source": [ "### Parameter der Simulation ###\n", "N = 10 # Anzahl der Teilchen\n", "box = (0.0, 1, 0.0, 0.5) # Größe der Box\n", "radius = 0.005 # Radius der Teilchen\n", "v = 0.4 # Anfangsgeschwindigkeit\n", "\n", "### Anfangsorte- und Geschwindigkeiten ###\n", "x, y = zeros(N), zeros(N)\n", "vx, vy = zeros(N), zeros(N)\n", "# Massen und Radien\n", "m, r = ones(N), radius * ones(N)\n", "initialize_particles_without_overlap!(x,y,vx,vy,r,v,box)\n", "\n", "Nmeasure = 10_000 # Anzahl der Messungen\n", "twait = 200 # Anzahl der Updateschritte zwischen den Messungen (groß für wirklich unabhängige Messungen)\n", "\n", "ps = [0.1,0.3,0.5,0.7,0.9] # Unterschiedliche Werte für p\n", "pleft = similar(ps) # Array um Anzahl der \"Teilchen links\"-Events zu speichern.\n", "\n", "# Schleife über alle p\n", "for i in eachindex(ps)\n", " left = 0 # Variable die Anzahl der Events \"Teilchen links\" zaehlt.\n", " # Schleife über alle Messungen\n", " for _ in 1:Nmeasure\n", " # Warte twait update Schritte zwischen zwei Messungen.\n", " for _ in 1:twait\n", " # Teilchen kollidieren.\n", " collide_particles!(x,y,vx,vy,m,r)\n", " # Teilchen fortbewegen.\n", " move_particles!(x,y,vx,vy,0.01,r,box)\n", " end\n", " \n", " # Messe am 1. Teilchen: Wenn im linken Teil mit Volumen ps[i]*box[2], erhöhe left\n", " (x[1] < ps[i] * box[2]) && (left += 1)\n", " end\n", " \n", " pleft[i] = left\n", "end \n" ] }, { "cell_type": "code", "execution_count": null, "id": "551dd900", "metadata": {}, "outputs": [], "source": [ "# Darstellung der Haeufigkeiten als barplot.\n", "fig, ax, p = barplot(ps, pleft)\n", "\n", "# Plot verschönern\n", "ax.xlabel = \"p\"\n", "ax.ylabel = \"# Teilchen links\"\n", "ax.yticks = 1000:2000:9000\n", "ax.xticks = 0.1:0.2:0.9\n", "\n", "fig" ] }, { "cell_type": "markdown", "id": "6e010d85", "metadata": {}, "source": [ "### Jetzt sind Sie gefragt: ###\n", "Genau wie oben wird jetzt wiederholt der Ort eines Teilchens gemessen. Zwischen jeder Messung wird erneut eine gewisse Zeit abgewartet. Diesmal führen wir aber viele *Messreihen* durch (z.B. 4000 Reihen), wobei jede Reihe aus einer kleineren Anzahl an Einzelmessungen (z.B. 25 Messungen) besteht. Für jede einzelne Messreihe sollen sie abspeichern, wie oft das Teilchen im linken Teil der Box mit Volumen $pV$ gemessen wird.\n", "Sie erhalten so 4000 Werte zwischen $k=0$ und $k=25$. Stellen Sie die Häufigkeiten der Werte von $k$ in einem Histogramm dar (mittels `hist`). Welcher Verteilung ähnelt das Histogramm? Stellen Sie diese Verteilung ebenfalls dar. " ] }, { "cell_type": "code", "execution_count": null, "id": "1e0e1d1b", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.8.5", "language": "julia", "name": "julia-1.8" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.8.5" } }, "nbformat": 4, "nbformat_minor": 5 }