Mathematisches Pendel animieren

[Marvin75854]

Hi,

ich habe ein kleines Programm geschrieben, welches ein einfaches Pendel mit verschiedenen numerischen Verfahren berechnet. Dieses Pendel will ich jetzt noch mithilfe von Python simulieren. Ich arbeite unter Python 2.7. Ich habe die Berechnung bereits durchgeführt und mit matplot geplottet. Habe jetzt schon eine menge gegooglet und was ausprobiert, aber irgendwie funktioniert nichts und ich verstehe nicht so wirklich wie das mit dem simulieren funktioniert. Der entsprechende Teil meines codes sieht so aus:

Code: Alles auswählen

import numpy as np
import matplotlib.pylab as plt
import math

...


class problem_2(object):
    def __init__(self, dt, frict):
        self.frict = frict 
        self.dt = dt
        self.theta0 = math.radians(85)
        self.dtheta0 = math.radians(0)
        self.m = 1.0
        self.g = 9.81
        self.l = 1.0
        self.t = 10.0
        self.end = int(self.t/self.dt)
        self.time = np.zeros(self.end + 1)
        self.theta = np.zeros(self.end + 1)
        self.dtheta = np.zeros(self.end + 1)
        self.time[0] = 0
        self.theta[0] = self.theta0
        self.dtheta[0] = self.dtheta0


    def computation_euler(self):
        for n in range (0, self.end):
            self.theta[n+1] = self.theta[n] + self.dt*self.dtheta[n]
            self.dtheta[n+1] = self.dtheta[n] + self.dt * ((-self.frict /self.m)*self.dtheta[n]+(-self.g/self.l) * math.sin(self.theta[n]))
            self.time[n+1] = self.time[n] + self.dt
        return (self.time, self.theta, self.dtheta)

    def plot(self, a):
        plt.figure()
        plt.plot(a[0], a[1]*180/math.pi)
        plt.xlabel("time [s]")
        plt.ylabel("angle [deg]")
        plt.grid()
        plt.show()


...
if __name__ == '__main__':
...
    y = problem_2
    euler_ts_1 = y.computation_euler(problem_2(0.001, 1.0))
    y.plot(problem_2(0.001, 1.0), euler_ts_1)

Man könnte die ganze Berechnung, bzw. den ganzen Code sicher auch besser aufbauen, aber das ist erstmal nicht so wichtig. Also das ganze will ich jetzt noch simulieren, habe aber absolut keinen Ansatzpunkt. Hoffe jemand kann mir helfen.

Gruß
Marvin

[__blackjack__]

@Marvin75854: Es gibt schon einen Unterschied zwischen „besser aufbauen“ und es nicht so völlig Banane zu schreiben. Die Klasse macht überhaupt keinen Sinn und wird auch komplett falsch verwendet. Das *ist* wichtig! Du hast Klassen sehr offensichtlich so überhaupt nicht verstanden. Warum sind das nicht einfach zwei Funktionen?

[__blackjack__]

Hier ist der Code ohne Klasse:

Code: Alles auswählen

import math
import matplotlib.pylab as plt
import numpy as np

#...


def computation_euler(dt, frict):
    theta0 = math.radians(85)
    dtheta0 = math.radians(0)
    m = 1.0
    g = 9.81
    l = 1.0
    t = 10.0
    end = int(t/dt)
    time = np.zeros(end + 1)
    theta = np.zeros(end + 1)
    dtheta = np.zeros(end + 1)
    time[0] = 0
    theta[0] = theta0
    dtheta[0] = dtheta0
    for n in range(end):
        theta[n+1] = theta[n] + dt*dtheta[n]
        dtheta[n+1] = dtheta[n] + dt * ((-frict /m)*dtheta[n]+(-g/l) * math.sin(theta[n]))
        time[n+1] = time[n] + dt
    return (time, theta, dtheta)


def plot(a):
    plt.figure()
    plt.plot(a[0], a[1]*180/math.pi)
    plt.xlabel("time [s]")
    plt.ylabel("angle [deg]")
    plt.grid()
    plt.show()


#...
if __name__ == '__main__':
#...
    plot(computation_euler(0.001, 1.0))

Als nächstes solltest Du dann teilweise schlechte Namensgebung in Angriff nehmen und Numpy dann auch mal tatsächlich verwenden statt da eine ``for``-Schleife selbst zu schreiben. Oder falls das nicht geht, einfach Listen verwenden und Werte dort anhängen statt Arrays mit 0en vor zu belegen und dann Index für Index durch den tatsächlichen Wert zu ersetzen. Das ist „unpythonisch“.

[Marvin75854]

@Marvin75854: Es gibt schon einen Unterschied zwischen „besser aufbauen“ und es nicht so völlig Banane zu schreiben. Die Klasse macht überhaupt keinen Sinn und wird auch komplett falsch verwendet. Das *ist* wichtig! Du hast Klassen sehr offensichtlich so überhaupt nicht verstanden. Warum sind das nicht einfach zwei Funktionen?

Ja so zu 100% hab ich Klassen auch nicht verstanden. In meinem Code zeige ich ja nur eine Berechnungsmethode, das Euler-Verfahren. Mein komplettes Programm, bzw. die Klasse beinhaltet aber noch andere Lösungsverfahren. Wollte jetzt nicht alles posten, weil es für meine eigentliche Frage nach einer Animation egal ist. Ich habe das also einfach in eine Klasse geschrieben um die ganzen Variablen nicht 5 mal definieren zu müssen. Ich hätte das auch global oder in der main-methode machen können, nehme ich an. Ohne Klasse wäre mir das schreiben des Programms vermutlich sogar deutlich leichter gefallen, dachte aber das wäre ganz gut, wenn ich das so machen würde.

Als nächstes solltest Du dann teilweise schlechte Namensgebung in Angriff nehmen und Numpy dann auch mal tatsächlich verwenden statt da eine ``for``-Schleife selbst zu schreiben. Oder falls das nicht geht, einfach Listen verwenden und Werte dort anhängen statt Arrays mit 0en vor zu belegen und dann Index für Index durch den tatsächlichen Wert zu ersetzen. Das ist „unpythonisch“.

Ok, vielleicht habe ich nicht ganz verstanden wie das mit numpy funktioniert. Ich habe einen Teil des Codes aus einer Vorlesung abgekupfert und für mich angepasst. Ich erzeuge doch erst einen Array mit der Länge end + 1, z.b. time = time und dann gebe ich doch nur dem ersten Wert des Arrays eine 0 oder? In der for-schleife nähere ich mich dann mit dem Euler-Verfahren einer Differenzialgleichung an. Ich berechne also basierend auf dem jetzigen Wert für z.b. time [n] den wert time [n+1]. und dann auf dem Wert für time [n+1] berechne ich den nächsten Wert time [n+2]. Wenn es dafür eine bessere Funktion als eine for-Schleife gibt, kenne ich diese nicht. Ich nehme an, dass wir das mit Arrays machen sollen und nicht mit Listen, da das ganze eben auch relativ schnell in Matrizenrechnung ausarten kann und vermutlich auch wird.
Die Namen sind vielleicht etwas verwirrend, ich habe aber einfach die Formelzeichen der Differentialgleichung verwendet. Also beispielsweise m für Masse. Ich persönlich komme damit gut klar.

EDIT: Ah durch time = np.zeros(end +1) fülle ich den gesamten Array mit 0en. Ich könnte das vermutlich auch einfach so machen time = np.array(end+1)? Habe den Teil ehrlich gesagt einfach so übernommen und mir da nicht so groß Gedanken drüber gemacht.

[__blackjack__]

@Marvin75854: Die Erklärung zum Sinn der Klasse verstehe ich nicht. Die wird auch bei mehr Verfahren nicht sinnvoller, also da braucht man IMHO nicht mehr vom Rest des Codes sehen, das ist einfach keine sinnvolle Klasse gewesen. Ich habe ja die beiden Funktionen gezeigt, die das letztendlich sind. `plot()` war von vornherein einfach nur eine Funktion auch wenn sie in der Klasse als ”Methode” stand und wenn man eine Klasse hat die nur aus einer `__init__()` und *einer* weiteren Methode besteht, hat man in der Regel eine Funktion unnötig kompliziert geschrieben. Ausnahme währen hier Klassen die hauptsächlich dazu dienen Werte zu einem Datensatz zusammen zu fassen, aber dieser Fall liegt hier nicht vor.

Numpy-Arrays verwendet man wenn man die speziellen Operationen auf diesen Arrays auch verwendet. Also wenn man keine Schleifen in Python schreibt, sondern Operationen und Methoden von Numpy verwendet, wo die Schleifen in dem Code dieser Array-Datentypen ablaufen.

`time` zum Beispiel würde man nicht so umständlich per Hand erstellen, und dazu auch noch mit wiederholten Additionen wo sich die Fehler durch die Gleitkommaungenaugkeiten aufsummieren, sondern einfach `numpy.linspace()` aufrufen.

[Marvin75854]

@Marvin75854: Die Erklärung zum Sinn der Klasse verstehe ich nicht. Die wird auch bei mehr Verfahren nicht sinnvoller, also da braucht man IMHO nicht mehr vom Rest des Codes sehen, das ist einfach keine sinnvolle Klasse gewesen. Ich habe ja die beiden Funktionen gezeigt, die das letztendlich sind. `plot()` war von vornherein einfach nur eine Funktion auch wenn sie in der Klasse als ”Methode” stand und wenn man eine Klasse hat die nur aus einer `__init__()` und *einer* weiteren Methode besteht, hat man in der Regel eine Funktion unnötig kompliziert geschrieben. Ausnahme währen hier Klassen die hauptsächlich dazu dienen Werte zu einem Datensatz zusammen zu fassen, aber dieser Fall liegt hier nicht vor.

Naja ich dachte nur, dass wenn ich 5 Methoden schreibe, die alle die gleichen Werte brauchen, würde es Sinn machen die alle zusammen in eine Klasse zu schreiben. Die Plot Funktion habe ich halt einfach noch mit reingeschrieben.

Numpy-Arrays verwendet man wenn man die speziellen Operationen auf diesen Arrays auch verwendet. Also wenn man keine Schleifen in Python schreibt, sondern Operationen und Methoden von Numpy verwendet, wo die Schleifen in dem Code dieser Array-Datentypen ablaufen.

`time` zum Beispiel würde man nicht so umständlich per Hand erstellen, und dazu auch noch mit wiederholten Additionen wo sich die Fehler durch die Gleitkommaungenaugkeiten aufsummieren, sondern einfach `numpy.linspace()` aufrufen.

Ok, time ist ein schlechtes Beispiel, weil da der Abstand zwischen den Punkten immer gleich ist. Bei den anderen beiden Arrays, muss der jeweils nächste Wert aber ja immer mit einer Differentialgleichung errechnet werden. Also um theta [n+1] zu berechnen, muss ich zuvor erstmal theta [n] berechnet haben. Und der Abstand ist dann eben jedes mal anders. Ich kenne wie gesagt nicht alle Funktionen von numpy, aber würde mich extrem wundern, wenn diese ganzen verschiedenen numerischen solver da schon zur Verfügung stehen. Ich denke auch, dass ich Numpy verwenden soll, einfach aus dem Grund, dass ich momentan nur eine Differentialgleichung betrachte, in zwei, drei Wochen habe ich jedoch ein ganzes System an Differentialgleichungen. Und dafür brauche ich dann diese ganzen Numpy funktionen. Also irgendwelche Matritzenmultiplikationen oder lineare Gleichungssysteme lösen usw.

Wie bekomme ich denn jetzt aus dem Ergebnis was ich habe eine Animation? Ich habe in dem Array theta[] jetzt eine Schwingung, also den Winkel des Pendels in Abhängigkeit von der Zeit. Damit müsste es doch irgendwie möglich sein das ganze zu animieren.

[__deets__]

Hast du mal matplotlib Animationen angeschaut? Die haben sogar ein Pendel-Beispiel...

https://matplotlib.org/examples/animati ... mated.html

[__blackjack__]

@Marvin75854: Wenn Du 5 Funktionen hast die die gleichen Argumente bekommen sollen und das relativ viele sind, macht es eher Sinn zu schauen ob man welche von diesen Argumenten sinnvoll zu einer (Daten)Klasse zusammenfassen kann. Wenn sie als Einheit unveränderlich sind und die Klasse sonst keine weiteren Methoden oder Properties bekommt, dann zum Beispiel als `collections.namedtuple()`.

Klassen sind nicht dazu da um es sich zu ersparen einen beliebigen Satz an Werten nicht immer übergeben zu müssen. Man muss der Klasse einen sinnvollen Namen geben können. Wenn man das nicht kann, hat man in der Regel Sachen verrührt die nicht zusammen gehören.

Für `theta` und `dtheta` sehe ich jetzt auch nicht auf anhieb wie das `numpy` Sinn machen würde. So eine Reihe wo die Werte immer auf den vorherigen Werten basieren würde ich als Generatorfunktion umsetzen und mir dann die entsprechende Anzahl von Elementen daraus nehmen, und dann beispielsweise in ein Numpy-Array stecken wenn das für die weitere Verarbeitung Sinn macht.

Die Rechnung ``* 180 / math.pi`` ist das Gegenteil von `radians()`, also `degrees()` und das kennt `numpy` auch für Arrays statt für Skalarwerte.

Man kännte dann beispielsweise ungefähr hier heraus kommen:

Code: Alles auswählen

#!/usr/bin/env python3
from collections import namedtuple
import math

import matplotlib.pylab as plt
from more_itertools import take
import numpy as np

GRAVITY = 9.81
#...

Pendulum = namedtuple('Pendulum', 'mass length start_theta start_dtheta')


def iter_euler(pendulum, dt, friction):
    theta = pendulum.start_theta
    dtheta = pendulum.start_dtheta
    while True:
        yield theta, dtheta
        theta, dtheta = (
            theta + dt * dtheta,
            dtheta
            + dt * ((-friction / pendulum.mass) * dtheta
            + (-GRAVITY / pendulum.length) * math.sin(theta))
        )


def computation_euler(pendulum, t, dt, friction):
    count = int(t / dt) + 1
    theta, dtheta = np.array(take(count, iter_euler(pendulum, dt, friction))).T
    return (np.linspace(0, t, count), theta, dtheta)


def plot(time, angle):
    plt.figure()
    plt.plot(time, np.degrees(angle))
    plt.xlabel('time [s]')
    plt.ylabel('angle [deg]')
    plt.grid()
    plt.show()


def main():
    #...
    pendulum = Pendulum(1, 1, math.radians(85), math.radians(0))
    time, theta, _ = computation_euler(pendulum, 10, 0.001, 1)
    plot(time, theta)
    #...


if __name__ == '__main__':
    main()