{ "cells": [ { "cell_type": "markdown", "id": "99836a8a", "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": "1b70a269", "metadata": {}, "source": [ "\n", "\n", "
Aufgabe 12: Satz von Liouville [8 + 2 Punkte, davon 4 Punkte für Teil b)]
" ] }, { "cell_type": "markdown", "id": "0a56c073", "metadata": {}, "source": [ "### b) ### \n", "Mit Hilfe dieses Notebooks können Sie die Phasenraumdynamik unterschiedlicher physikalischer Systeme untersuchen. Genauer wollen wir die zeitliche Entwicklung eines *Phasenraumvolumens* nachvollziehen. In den ersten Zellen werden zunächst einige Funktionen bereit gestellt, die die gewünschten Anfangsbedingungen erzeugen und die die zu lösende Differentialgleichungssysteme integrieren. Sie können diese Funktionen unverändert übernehmen. Lesen Sie sich die nachfolgenden Zellen trotzdem aufmerksam durch, und ergänzen Sie anschließend weiter unten, an den markierten Stellen, Ihren Code, um die auf dem Aufgabenblatt gegebenen Systeme zu implementieren." ] }, { "cell_type": "markdown", "id": "015e7fc1", "metadata": {}, "source": [ "#### Einbinden der benötigten Pakete ####" ] }, { "cell_type": "code", "execution_count": null, "id": "6997cd3b", "metadata": {}, "outputs": [], "source": [ "using DifferentialEquations, GLMakie" ] }, { "cell_type": "markdown", "id": "a4a7704a", "metadata": {}, "source": [ "#### Integration der Differentialgleichung um einen Zeitschritt. ####" ] }, { "cell_type": "code", "execution_count": null, "id": "1ce4d74a", "metadata": {}, "outputs": [], "source": [ "# Funktion integriert die Differentialgleichung um den Zeitschritt dt weiter.\n", "# Hierbei ist \n", "# - integ ein ODEIntegrator (spezieller Datentyp aus dem Paket DifferentialEquations)\n", "# - points eine Observable, die die aktuellen Orte und Impulse speichert.\n", "# - dt der Zeitschritt, den die DGL nach vorne entwickelt wird.\n", "function animstep!(integ, points, dt)\n", " \n", " # Differentialgleichung einen Schritt weiter integrieren.\n", " step!(integ, dt, true)\n", " \n", " # Update der Observablen, die die aktuellen Orte/Impulse speichert.\n", " for i in 1:size(points[], 1)\n", " points[][i] = Point2f(integ[2*i-1], integ[2*i])\n", " end\n", " \n", " # Observable auf Update aufmerksam machen\n", " notify(points)\n", " \n", "end" ] }, { "cell_type": "markdown", "id": "a62ea31a", "metadata": {}, "source": [ "#### Erzeuge Anfangsbedingungen. ####" ] }, { "cell_type": "code", "execution_count": null, "id": "b4042846", "metadata": {}, "outputs": [], "source": [ "# Erzeuge zufällige Anfangsbedingungen (x,p), die gleichmäßig verteilt innerhalb eines Kreises \n", "# mit Mittelpunkt (x0,p0) und Radius R liegen. N ist die Anzahl der zurückgegebenen Wertepaare.\n", "function get_init_inner(x0,p0,R,N)\n", "\n", " # generate equally distributed random initial condition in circle.\n", " r, phi = R * sqrt.(rand(N)), rand(N) * 2 * pi\n", " x, p = x0 .+ r .* cos.(phi), p0 .+ r .* sin.(phi)\n", "\n", " # provide correct formatting: [q1,p1,...,qN,pN]\n", " u0 = zeros(2*N)\n", " for i in 1:N\n", " u0[2*i-1] = x[i]\n", " u0[2*i] = p[i]\n", " end \n", " \n", " return u0\n", "end\n", "\n", "# Erzeuge Anfangsbedingungen (x,p), die gleichmäßig verteilt auf einem Kreis mit Mittel-\n", "# punkt (x0,p0) und Radius R liegen. N ist die Anzahl der zurückgegebenen Wertepaare\n", "function get_init_border(x0,p0,R,N)\n", "\n", " # generate equally distributed initial condition on circle.\n", " phi = range(0,2*pi, length = N)\n", " x, p = x0 .+ R .* cos.(phi), p0 .+ R .* sin.(phi)\n", "\n", " # provide correct formatting: [q1,p1,...,qN,pN]\n", " u0 = zeros(2*N)\n", " for i in 1:N\n", " u0[2*i-1] = x[i]\n", " u0[2*i] = p[i]\n", " end \n", " \n", " return u0\n", "end" ] }, { "cell_type": "markdown", "id": "cee58e21", "metadata": {}, "source": [ "#### Definiere das DGL Problem. ####" ] }, { "cell_type": "code", "execution_count": null, "id": "cc4f7d6b", "metadata": {}, "outputs": [], "source": [ "# Funktion erzeugt zwei ODE Probleme (für das Innere und den Rand des Volumenelements)\n", "# Ninner / Nborder sind die jeweilige Anzahl der Anfangspunkte, R ist der Kreisradius.\n", "# Die Funktion problem definiert das Differentialgleichungsproblem (siehe unten).\n", "function generate_ODEprob(Ninner, Nborder, R, problem::Function)\n", " \n", " tspan = (0,10)\n", " params1, params2 = [Ninner], [Nborder]\n", " \n", " u0inner = get_init_inner(0,0,R,Ninner)\n", " u0border = get_init_border(0,0,R,Nborder)\n", "\n", " prob_inner = ODEProblem(problem, u0inner, tspan, params1)\n", " prob_border = ODEProblem(problem, u0border, tspan, params2)\n", "\n", " # initialize integrators for two sets of solutions\n", " integ_inner = init(prob_inner, save_everystep = false);\n", " integ_border = init(prob_border, save_everystep = false);\n", "\n", " # Observablen zur Darstellung.\n", " points_inner = [Point2f(u0inner[2*i-1], u0inner[2*i]) for i in 1:Ninner]\n", " points_border = [Point2f(u0border[2*i-1], u0border[2*i]) for i in 1:Nborder]\n", " points_inner = Observable(points_inner)\n", " points_border = Observable(points_border)\n", " \n", " return integ_inner, integ_border, points_inner, points_border\n", "end " ] }, { "cell_type": "markdown", "id": "0893d392", "metadata": {}, "source": [ "### In den folgenden Zellen sind Sie gefragt: ###\n", "Lesen Sie sich die nachfolgende Zelle gründlich durch und vollziehen Sie nach, wie die Bewegungsgleichungen des Systems aus $N$ ungekoppelten, harmonischen Oszillatoren implementiert wurde. Führen Sie anschließend die übernächste Zelle aus und machen Sie sich mit der Simulation vertraut. Sie können per Mausklick ein neues Phasenraumvolument auswählen und die Zeitentwicklung erneut starten. Zoomen ist per Mausrad (oder äquivalenter Funktion) möglich, der sichtbare Bereich lässt sich mit gedrückter rechter Maustaste verschieben.\n", "\n", "Ergänzen Sie anschließend - analog zum harmonischen Oszillator - die Differentialgleichungen, die das gedämpfte Pendel und das System mit der auf dem Aufgabenblatt gegebenen Hamiltonfunktion beschreiben. Für das gedämpfte Pendel können Sie $\\gamma = 0.1$, $m =1$ und $ K = 3$ setzen.\n", "Visualisieren Sie anschließend auch die Phasenraumdynamik der beiden neuen Systeme und beschreiben Sie Ihre Beobachtungen. \n", "\n", "**Optional:** Implementieren Sie auch die Bewegungsgleichungen des Duffing-Oszillators auf geeignete Weise und untersuchen Sie das System für verschiedene Parametersätze, etwa $\\alpha = -1$, $\\beta = 1$, $\\gamma = 0.02$, $\\delta = 3$, $\\omega = 1$." ] }, { "cell_type": "code", "execution_count": null, "id": "6cf7432b", "metadata": {}, "outputs": [], "source": [ "# Beispiel fuer Implementierung der DGL: \n", "# Ein einzelner Oszillator.\n", "# du = [dq/dt, dp/dt]\n", "# u = [q,p]\n", "# p = [sonstige Parameter]\n", "# t = Zeit (falls DGL explizit zeitabhängig wäre)\n", "function oscillator(du, u, p, t)\n", " du[1] = u[2]\n", " du[2] = -u[1]\n", "end\n", "\n", "# N ungekoppelte Oszillatoren.\n", "# du = [q1,p1,...,qN,pN]\n", "# du = [dq1/dt, dp1/dt, ... , dqN/dt, dpN/dt]\n", "# p = [sonstige Parameter]. Hier: Anzahl der Oszillatoren.\n", "# t = Zeit (falls DGL explizit zeitabhängig wäre)\n", "function Noscillators(du, u, p, t)\n", " N = p[1]\n", " for i in 1:2:(2*N-1)\n", " du[i] = u[i+1]\n", " du[i+1] = -u[i]\n", " end\n", " return\n", "end\n", "\n", "\n", "function Ndamped_oscillators(du, u, p, t)\n", " ### Ergänzen Sie hier Ihren Code ###\n", " return\n", "end\n", "\n", "### Ergänzen Sie hier die Bewegungsgleichungen der gegebenen Hamilton-Funktion ###" ] }, { "cell_type": "code", "execution_count": null, "id": "819c1ba0", "metadata": {}, "outputs": [], "source": [ "# Definiere das DGL Problem:\n", "\n", "# Anzahl an Punkten innen und auf Rand, Radius der Anfangsbedingungen.\n", "NI, NB, R = 1000, 1001, 5\n", "\n", "# Initialisiere das DGL Problem und die Observablen zur Darstellung der Trajektorien.\n", "IInner, IBorder, PInner, PBorder = generate_ODEprob(NI, NB, R, Noscillators)\n", "\n", "\n", "########################################################\n", "######### Visualisierung mit Makie #####################\n", "#########(kann unverändert übernommen werden) ########## \n", "########################################################\n", "\n", "fig = Figure()\n", "ax = Axis(fig[1,1])\n", "scatter!(ax, PInner)#, color = xvals, colormap = :viridis)\n", "scatter!(ax, PBorder)#, color = xvals, colormap = :viridis)\n", "\n", "# Makie Einstellungen.\n", "limits!(ax, -8,8,-8,8)\n", "hidedecorations!(ax)\n", "Makie.deactivate_interaction!(ax, :rectanglezoom)\n", "spoint = select_point(ax.scene, marker = :circle)\n", "\n", "display(fig)\n", "\n", "# Start/Stop Button\n", "run = Button(fig[2,1]; label = \"Start / Stop\", tellwidth = false)\n", "isrunning = Observable(false)\n", "on(run.clicks) do clicks; isrunning[] = !isrunning[]; end\n", "on(run.clicks) do clicks\n", " @async while isrunning[]\n", " isopen(fig.scene) || break\n", " \n", " # Zeitentwicklung beider DGLs.\n", " animstep!(IInner, PInner, 0.025)\n", " animstep!(IBorder, PBorder, 0.025)\n", " \n", " sleep(0.01)\n", " end\n", "end \n", "\n", "# Wähle neue Anfangsbedingungen per Mausklick.\n", "on(spoint) do z \n", " x,y = z\n", " \n", " # Neue Anfangsbedingungen\n", " uin = get_init_inner(x,y,R,NI)\n", " ubor = get_init_border(x,y,R,NB)\n", "\n", " # Übergabe der Anfangsbedingungen an Integrator.\n", " reinit!(IInner, uin)\n", " reinit!(IBorder, ubor)\n", " \n", " # Observable Updaten.\n", " PInner[][:] = [Point2f(IInner[2*i-1], IInner[2*i]) for i in 1:NI]\n", " PBorder[][:] = [Point2f(IBorder[2*i-1], IBorder[2*i]) for i in 1:NB]\n", " \n", " # Observable auf Update aufmerksam machen\n", " notify(PInner)\n", " notify(PBorder)\n", "end\n" ] }, { "cell_type": "code", "execution_count": null, "id": "f2e9a936", "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 }