Simulation von Differentialgleichungen mit Tabellenkalkulationen (2)

In diesem Beitrag soll es nun um das sogenannte Fadenpendel (oder auch mathematische Pendel) gehen – und wie man es mit Hilfe von OpenOffice „lösen“ kann! Dazu zunächst ein wenig Theorie zum Fadenpendel, bevor dann eine verblüffend einfache Lösung mittels Tabellenkalkulation folgt. Notwendig zum Verständnis ist auf jeden Fall mein vorheriger Artikel, in dem ich das grundlegende Prinzip der Näherungslösungen beschreibe. Hier möchte ich nur ganz kurz hinzufügen, wie mit einer zweiten Ableitung umgegangen wird.

Die Zweite Ableitung

Die 2. Ableitung einer Funktion ist – wer hätte das gedacht – die Ableitung der 1. Ableitung einer Funktion. Richtig schreiben müsste man für die zweite Zeitableitung also:

[latex size=“1″]ddot{f} = frac{d dot{f}}{dt}[/latex]

Wendet man nun den Gedanken aus dem vorherigen Artikel auch hierauf an, so kann ich – für hinreichend kleine [latex size=“1″]Delta t[/latex] die Ableitung näherungsweise berechnen über

[latex size=“1″]ddot{f} = frac{Delta dot{f}}{Delta t}[/latex]

Die Änderung der 1. Ableitung [latex size=“1″]Delta dot{f}[/latex] lässt sich also direkt durch Umstellen der Gleichung berechnen. Neu bei der Behandlung von Differentialgleichungen ist hier nur, dass auch für die 1. Ableitung ein Anfangswert vorgegeben werden muss. Dieser kann aber auch als 0 angenommen werden. Dies ist bereits alles an „neuem“ Hanwerkszeug, was wir für das Fadenpendel benötigen.

Das Fadenpendel

Jeder von uns kennt es: Ein schwerer Gegenstand, zum Beispiel eine Uhr, ist an einem Faden (oder einer Kette) angehängt. Wird der Gegenstand leicht ausgelenkt, so beginnt das Pendel hin- und her zu schwingen. Oder kurz gesagt: Es pendelt.

Um zur Differentialgleichung des Fadenpendels zu gelangen, müssen wir eine kleine Vorüberlegung anstellen. Am einfachsten ist die Kräftebetrachtung, wenn das Problem in ebenen Polarkoordinaten formuliert wird. Das hört sich schlimmer an, als es ist: Als Paramter verwende ich nicht x und y (das wäre ein karthesisches Koordinatensystem), sondern den Winkel [latex size=“1″]phi (t)[/latex], um den das Pendel zu einer bestimmten Zeit t ausgelenkt ist. Dieser hängt mit dem Weg (in karthesischen Koordinaten) über folgende Beziehung zusammen:

[latex size=“1″] phi (t) = s(t)/l[/latex]

Plausibel wird diese Beziehung bei der Betrachtung des folgenden Bildes aus der Wikipedia:

Der zu einem bestimmten Zeitpunkt zurückgelegte Weg s(t) ist – übertragen auf das Bild – gerade b, die Länge unseres Pendels r und der Auslenkungswinkel [latex size=“1″]phi[/latex] ist gerade [latex size=“1″]alpha[/latex]. Hieraus folgt auch, dass die Beschleunigung a – die ja selbst die zweite Ableitung von s(t) nach der Zeit ist – proportional zur 2. Ableitung der Winkelfunktion nach der Zeit ist:

[latex size=“1″]ddot{s}(t)=l cdot ddot{alpha}[/latex].

Was hilft das nun? Mit Hilfe dieser Beziehung lässt sich die Grundgleichung der Mechanik

[latex size=“1″]F=m cdot a = m cdot ddot{s}[/latex]

umschreiben in

[latex size=“1″]F=m cdot l cdot ddot{phi}[/latex].

In Worten bedeutet dies, dass die Kraft F gerade proportional zur Beschleunigung a ist (mit dem Proportionalitätsfaktor „Masse“ m) – und auch proportional zur zweiten Ableitung des jeweiligen Auslenkungswinkels [latex]phi[/latex] nach der Zeit, da die Beschleunigung gerade die zweite Ableitung der Streckenfunktion nach der Zeit ist. Gewonnen haben wir bis hierher eigentlich noch nicht viel.

Einen etwas intensiveren Blick müssen wir nun auf die beim Fadenpendel auftretenden Kräfte riskieren. Plausibel werden die folgenden Formeln hoffentlich mit folgender formschönen Skizze:

Die Komponentenzerlegung der Gewichtskraft am Fadenpendel

Dies ist nur die Gewichtskraft mg senkrecht nach unten. Diese lässt sich nun in eine Komponente in Fadenrichtung aufteilen ([latex size=“1″]F_{r}[/latex]) und in eine dazu senkrechte ([latex size=“1″]F_{T}[/latex]). Erstere ist genauso groß wie die Kraft, die der Faden auf die Masse ausübt und sorgt dafür, dass der Faden gespannt bleibt. Letztere ist die Kraft, die die Masse zum Schwingen bringt. In der Zeichnung erkennt man folgenden Zusammenhang zwischen Auslenkungswinkel [latex size=“1″]phi[/latex], der Beschleunigung a, der Erdbeschleunigung g und der Masse m:

[latex size=“1″]sin phi = frac{F_{T}}{m cdot g}[/latex] -> [latex size=“1″]F_{T} =- m cdot g cdot sin phi[/latex]

Das Minuszeichen rührt daher, dass die Kraft der Auslenkung entgegen wirkt. Nun muss noch die oben abgeänderte Grundgleichung mit diesem Ausdruck verheiratet werden und wir sind schon am Ende angelangt. Dieses erreichen wir durch Gleichsetzen:

[latex size=“1″]m cdot l cdot ddot{phi} =- m cdot g cdot sin phi[/latex]

Durch Kürzen der Masse und eine kleine Umstellung der Gleichung erhalten wir dadurch die Differentialgleichung des Fadenpendels. In dieser taucht die Masse m nicht mehr auf (die Schwingungsperiode ist unabhängig von dieser) und es handelt sich um eine lineare homogene Differentialgleichung 2. Ordnung (wenn der Luftwiderstand und Reibungsverluste vernachlässigt werden). Sie lautet:

[latex size=“1″]ddot{phi}=- frac{g}{l} cdot sin(phi)[/latex]

Für kleine Winkel (so im Bereich bis 10°) ist der Sinus des Winkels in sehr guter Näherung gleich dem Winkel selbst. Somit vereinfacht sich die Gleichung nochmal zu:

[latex size=“1″]ddot{phi}=- frac{g}{l} cdot phi[/latex]

Diese Differentialgleichung gilt es nun zu lösen – nach der im letzten Artikel beschriebenen Methode mit Hilfe von OpenOffice!

Lösung der Differentialgleichung mit OpenOffice

Die Orginaldatei mit meiner Simulation findet sich hier: fadenpendel. Ich werde mich im Folgenden auf die entsprechende Tabelle beziehen.

Zunächst wurden wieder die benötigten Werte als Felder definiert. Dieses sind die Anfangsauslenkung (Feld B4), die Erdbeschleunigung g (B5) die Länge des Pendels (B6) und das Diskretisierungsintervall (B7). Außerdem wird noch die Anfangsgeschwindigkeit (also der Anfangswert der ersten Ableitung von [latex size=“1″]phi[/latex]) benötigt (B8).

Als Tabellenspalten benötigt werden die Zeit (Spalte A), die Änderung der ersten Ableitung (Spalte B – in der Tabelle ist diese mit Deltaomega benannt, da in der Physik die Winkelgeschwindigkeit [latex size=“1″]omega = dot{phi}[/latex] eine gebräuchliche Abkürzung für die 1. Ableitung der Winkelfunktion nach der Zeit ist), den Wert der ersten Ableitung (Spalte C) und den Wert des Auslenkungswinkels (Spalte D).

Zunächst erhalten alle ihre Anfangswerte (Zeile 10). Die Zeit wird dann um das feste Intervall dt inkrementiert. Die Änderung der ersten Ableitung ergibt sich, wenn die zweite Ableitung geschrieben wird als

[latex size=“1″]frac{Delta omega}{Delta t} = – frac{g}{l} cdot phi[/latex]

durch Umstellen der obigen Gleichung. Dieses ist im Feld B11 implementiert, wobei davon ausgegangen wird, dass der Auslenkungswinkel in D10 über diesen (kurzen) Zeitraum nahezu konstant geblieben ist. Die Änderung, die in B11 berechnet wurde, wird nun in C11  mit dem vorherigen Wert C10 der 1. Ableitung verrechnet. In guter Näherung ergibt sich der neue Auslenkungswinkel in D11 dann durch die Summe aus vorheriger Auslenkung in D10 und dem Produkt von 1. Ableitung und eingestelltem Zeitintervall.

Und das war es auch schon. Bei mir entsteht dadurch dieser wunderhübsche, sinusförmige Graph. Gut zu sehen ist auch die Phasenverschiebung zwischen Winkelgeschwindigkeit omega und Auslenkungswinkel phi (nämlich gerade 90°).



Simulation des Fadenpendels mit OpenOffice
Simulation des Fadenpendels mit OpenOffice

Wozu das ganze?

Das Fadenpendel ist in der Oberstufenphysik ein häufig diskutiertes Modellsystem. Jedoch wird die Lösung der Differentialgleichung mehr als trigonometrische Funktion „geraten“ denn wirklich gerechnet. Dieses intuitive Raten hat seine Berechtigung und hilft natürlich auch bei der Diskussion des Harmonischen Oszillators (dessen DGL im reibungsfreien Fall strukturell genauso aussieht wie diejenige des Fadenpendels).

Wird jedoch mit einem Oberstufenkurs auch die numerische Lösung in Excel (zusätzlich) durchgeführt, ergeben sich in meinen Augen viele Vorteile:

  • die Schüler erwerben Kenntnisse im Umgang mit Excel – oder eben auch OpenOffice. Diese sind auch außerhalb des physikalischen Kontextes wichtig.
  • die Schüler konstruieren sich – natürlich nach Anleitung – mit vergleichsweise wenig Mathematik die Lösung einer auf den ersten Blick komplizierten Gleichung selbst.
  • die Aufgabe eigent sich wunderbar, um in Partnerarbeit durchgeführt zu werden. Nach einer kurzen Einführung durch den Lehrer (der beispielsweise die Simulation einer DGL 1. Ordnung vorführt) können die Schüler selbstständig in einer sozialen Lernsituation eine Lösung erarbeiten.
  • der Umgang mit dem Computer hat eine große motivationale Komponente.

Im nächsten Artikel will ich mit der hier vorgestellten Methode auch den harmonischen Oszillator noch näher unter die Lupe nehmen. Kommentare und Verbesserungsvorschläge sind jederzeit willkommen!

Bisher in dieser Artikelserie:

Simulation von Differentialgleichungen mit Tabellenkalkulationen (1)

Differentialgleichungen bestimmen die Physik: Von den Newton’schen Bewegungsgleichungen bis hin zur Schrödingergleichung der Quantenmechanik. Leider ist die zur Lösung benötigte Mathematik nicht immer ganz einfach. Und schon gar nicht immer anschaulich. Vor allem in der Schule stellen DGL’s Schüler oft vor Verständnisschwierigkeiten.

Moderne Computer bieten hier viele Möglichkeiten, den Differentialgleichungen ihren Schrecken zu nehmen und Differentialgleichungen – gerade für Schüler – anschaulich zu machen. Nummerische (Näherungs-)lösungen sind ohne großen Aufwand mit Standardsoftware möglich, oftmals trägt das selbstständige „Programmieren“ sogar erheblich zum Verständnis der Physik hinter den DGL bei. Ich möchte hier nach und nach einige Aspekte etwas näher beleuchten und Zugänge  mit Excel OpenOfficeCalc, vorstellen. Heute geht es erstmal um ein wenig Theorie, die Entladung eines Kondensators und Reaktionen 1. Ordnung als Anwendung für die chemische Kinetik.

Näherungslösungen für Differentialgleichungen

Eine Differentialgleichung stellt eine Beziehung zwischen einer Funktion und ihren Ableitungen her (im einfachsten Fall: zwischen Funktion und ihrer ersten Ableitung). Die Lösungsmenge so einer Differentialgleichung sind dann all diejenigen Funktionen, die die in der DGL gestellten Bedingungen erfüllen. Ein einfaches Beispiel (das uns heute in abgewandelter Form beschäftigen wird) ist eine lineare, homogene DGL 1. Ordnung:

[latex size=“1″]frac{dx}{dt}=x leftrightarrow dot{x} = x[/latex]

Ich benutze hierbei die in der Physik übliche Schreibweise mit einem Punkt über der Variable, um eine Zeitableitung anzudeuten. Diese DGL ist analytisch einfach lösbar, auch viele Schüler können leicht nachvollziehen, dass gerade die e-Funktion ihre eigene Ableitung ist. Dennoch kann man hieran einige einfache Überlegungen anschließen. Zunächst noch einmal zu den Begrifflichkeiten: linear heißt eine DGL, wenn die Funktion und ihre Ableitungen nur mit dem Exponenten 1 in der Gleichung auftauchen. homogen, wenn es keinen Term gibt, in dem weder die Funktion noch eine ihrer Ableitungen auftauchen (entschuldigt, das ist mathematisch alles nur halbsauber, muss uns hier aber nicht näher interessieren).

Wie diskretisiert man nun diese DGL? Die Ableitung ist ja eigentlich nichts anderes als der folgende Grenzwert (der – so die Funktion differenzierbar ist – gegen einen bestimmten Wert konvergiert):

[latex size=“1″]limlimits_{Delta t rightarrow 0}{frac{Delta x}{Delta t}}[/latex]

Für kleine [latex]Delta t[/latex] ist aber bereits der Bruch eine sehr gute Näherung für die Ableitung – und genau das ist der Kerngedanke bei der numerischen Simulation. Kurz: Wähle [latex]Delta t[/latex] hinreichend klein, dann ist der Differenzenquotient

[latex size=“1″]frac{Delta x}{Delta t} = frac{dx}{dt}[/latex]

Das Gleichheitszeichen ist mathematisch natürlich nicht ganz sauber, der Differenzenquotient ist nur im Grenzwert gleich dem Differentialquotienten. Für unsere Zwecke ist es aber hilfreich, sich dieses Gleichheitszeichen zu denken.

DGL des Entladevorgangs eines Kondensators

Zur Herleitung der Differentialgleichung für den Entladevorgang eines Kondensators bietet die Wikipedia wesentlich mehr Informationen, als ich hier unterbringen kann und möchte. Entscheidend ist: Der Entladevorgang eines Kondensators wird durch folgende Differentialgleichung beschrieben:

[latex size=“1″]frac{dQ}{dt} =(-1) cdot frac{Q}{RC} cdot Q(t)[/latex]

C und R sind hierbei die bauartbedingte Kapazität und der ohm’sche Widerstand des Kondensators.

Simulation des Entladevorgangs mit OpenOfficeCalc

Eine von mir erstellte Beispieldatei gibt es hier. Auf diese wird im weiteren Text immer wieder Bezug genommen.

Notwendig zur Lösung der DGL ist zunächst die Diskretisierung: Der Differentialquotient [latex]frac{dQ}{dt}[/latex] wird durch den Differenzenquotienten [latex]frac{Delta Q}{Delta t}[/latex] angenähert (siehe oben). Weiterhin muss ein Startwert vorgegeben werden, der die Ladung zum Startzeitpunkt beschreibt. Zweckmäßiger Weise wird sowohl der Startwert als auch das Diskretisierungsintervall [latex]Delta t[/latex] über eine eigene Tabellenzeile definiert (in der Beispieldatei: Zellen B4 und B7).

Außerdem erscheint es mir zweckmäßig, den Widerstand und die Kapazität des Kondensators in einer eigenen Zelle zu definieren (B5 und B6). Damit ist auch schon fast alles getan und OpenOffice kann für uns rechnen.

Dazu wird zunächst eine Spalte mit fortlaufender Zeit definiert. Hierzu wird auf den Startwert „0“ in jeder Zeile das Diskretisierungsintervall [latex]Delta t[/latex] addiert (Spalte A ab A10). Das war noch einfach, oder?

Dann gilt es, die Ladungsänderung in dieser Zeit zu berechnen. Hierzu wird angenommen, dass die Ladungsmenge Q in der Zeit [latex]Delta t[/latex] als konstant angenommen werden kann. Dann gilt für [latex]Delta Q[/latex]:

[latex size=“1″]Delta Q =(-1) cdot frac{Q}{RC} cdot Delta t[/latex]

Dieser Wert wird das erste Mal zur Zeit [latex]Delta t[/latex] berechnet (B11) und gibt gerade die Änderung der Ladung an. Diese wird nun einfach zur vorher vorhandenen Ladung addiert (C11). Und das war es auch bereits: Die Zeile (A/B/C 11) kann einfach nach unten gezogen werden und das Ergebnis in ein Diagramm gezeichnet werden. Das sieht dann so aus:

Entladung eines Kondensators, simuliert mit OpenOffice
Entladung eines Kondensators, simuliert mit OpenOffice

Der Vorteil daran, die einzelnen Variablen als eigenes Feld in der Tabelle zu definieren: Es ist nun möglich, an den Werten herumzuspielen und zu schauen, welche Auswirkungen dies hat. Übrigens: Der Graph sieht doch verdächtig nach einer Exponentialfunktion aus, oder? Genau diese erhielte man auch durch die analytische Lösung der Differentialgleichung.

Simulation einer Reaktion 1. Ordnung

Nicht nur in der Physik tauchen solch einfache Differentialgleichungen 1. Ordnung auf, sondern z.B. auch in der chemischen Kinetik. Eine Reaktion 1. Ordnung gehorcht der folgenden DGL (c: Konzentration des Stoffes, k: Geschwindigkeitskonstante):

[latex size=“1″]frac{dc}{dt}=-k cdot c[/latex]

Sie sieht also auch ganz ähnlich aus wie oben. Die Lösung funktioniert daher ebenfalls analog und ich will sie gar nicht näher erläutern. Eine Beispieldatei befindet sich hier. Die Simulation gilt genauso übrigens auch für radioaktive Zerfälle, da diese im Grunde auch „nur“ Reaktionen 1. Ordnung darstellen…

Im nächsten Teil dieser kleinen Serie möchte ich versuchen, auch etwas kompliziertere Differentialgleichungen mit der hier vorgestellten Methode zu lösen. Zunächst werde ich mich an das mathematische Fadenpendel herantrauen, um dann auch noch den harmonischen Oszillator (mit Dämpfung!) zu besprechen. Kommentare und Ergänzungen sind natürlich immer willkommen!

Bisher in dieser Artikelserie: