Doppelpendel Differentialgleichung lösen

pleasediaaa · Gast

Meine Frage:
Guten Tag :)

ich bin Physikstudentin im ersten Semester und wir habe in unserem Computerkurs die Hausarbeit bekommen in Python unter anderem die Differentialglichung für ein Doppelpendel zu implementieren und wahlweise mit dem Eulerverfahren oder dem 4 Schritt Runge Kutter Verfahren zu lösen. Das Problem liegt momentan darin, dass wir in Mathe noch keine Differentialgleichungen hatten und uns deshalb 3 Tage lang was angelesen haben, in diesem Punkt jedoch verzweifeln.

Meine Ideen:
Das sind unsere Ideen: Bisher konnten wir mithilfe einiger Beiträge die Bewegungsgleichungen für die beiden Massen des Doppelpendels herleiten und auch verstehen. An Hand eines einfachen Pendels konnten wir auch schon Teile unseres Programmcodes testen und wissen,dass diese funktionieren (Eulerstep, Flowmap). Wie der Name schon sagt soll Eulerstep die differentialgleichung nach dem Eulerverfahren lösen. Flowmap soll den Pointcareschnitt für die Bedingungen, die in der Funktion pointcare aufgelistet sind berechnen.
Die Funktion startcondition soll zufällige Anfangswerte für die Geschwindigkeiten generieren. Die gleichung in startcondition ist abhängig von der durch den User eingegeben Systemenergie und der Bedingung, dass x1 und x2 0 sind (also nicht ausgelenkt)
x entspricht unserem Winkel
y der ersten Ableitung, also der Winkelgeschwindigkeit
und z der Winkelbeschleunigung
die indizes 1 und 2 jeweils für die Masse 1 oder Masse 2.
Unser Problem ist jetzt, dass die Bewegungsgleichungen aus der Funktion Doppelpendel ja voneinander abhängen. Das heißt will ich z1 haben, brauche ich z2 und andersrum. So wie wir es bisher beurteilen können, liegt darin unser Problem. Also eher ein mathematisches. Nichts destotroz, schreibe ich euch unseren Programmcode mit hin.
Im folgenden unser Programmcode:

import numpy as np
from pylab import*
%matplotlib inline

def doppelpendel(z1, y1, y2, t, dt, lam, mu, omega, e, x1, x2):
z2 = -lam*omega**2*np.sin(x2)+lam*y1**2*np.sin(x1-x2)-lam*z1*np.cos(x1-x2)
z1 = (omega**2)*np.sin(x1) - (mu/lam)*y2**2*np.sin(x1-x2) - (mu/lam)*z2*np.cos(x1-x2)
return z1, z2
def doppelpendelz1(z2, y1, y2, t, dt, lam, mu, omega, e, x1, x2):
z1 = (omega**2)*np.sin(x1) - (mu/lam)*y2**2*np.sin(x1-x2) - (mu/lam)*z2*np.cos(x1-x2)
print("z1 ist {}".format(z1))
return z1

def doppelpendelz2(y1, y2, t, dt, lam, mu, omega, e, x1, x2):
z2 = (-lam*omega**2*np.sin(x2) + lam*y1**2*np.sin(x1-x2) - omega**2*lam*np.cos(x1-x2)*np.sin(x1) + mu*y2**2*np.sin(x1-x2)*cos(x1-x2))/(1-mu*(np.cos(x1-x2))**2)
print("z2 ist {}".format(z2))
return z2

def startcondition(lam, mu, e, steps):
cachey1 = []
cachey2 = []
for i in range(1, steps+1):
y1 = np.random.uniform(-100,100)
y21 = -lam*y1 + sqrt(lam**2 * y1**2 - (lam**2 * y1**2 - 2*lam*e)/ mu)
y22 = -lam*y1 - sqrt(lam**2 * y1**2 - (lam**2 * y1**2 - 2*lam*e)/ mu)
cachey1.append(y1)
cachey2.append(y21)
cachey1.append(y1)
cachey2.append(y22)
return y1, y21

def EulerStep(y1, y2, f1, f2, dt,t, lam, mu, omega, e, x1, x2):
z2 = f1(y1, y2, t, dt, lam, mu, omega, e, x1, x2)
z1 = f2(z2, y1, y2, t, dt, lam, mu, omega, e, x1, x2)

y2 += dt*z2
x2 += dt*y2

y1 += dt*z1
x1 += dt*y1
t += dt
return x1, x2, y1, y2, t, dt, z1, z2

def FlowMap2(dt,t, steps, f1, f2, Integrator, stopcond, lam, mu, omega, e):
cachey2=[]
cachex2=[]
cachet=[]
cachez1=[]
y1, y21 = startcondition(lam, mu, e, steps)
j=1
y2 = y21
x1=0
x2=0
print( "lam = {} mu = {} omega = {} e = {} y1 = {} y2 = {}".format(lam,mu,omega,e,y1,y2))
for i in range(1, steps+1):
cond =False
while cond==False:
x1, x2, y1, y2, t, dt,z1,z2 = Integrator(y1, y2, f1, f2, dt,t, lam, mu, omega, e, x1, x2)
cond = stopcond(x1,y1)
cachez1.append(z1)
cachet.append(t)
cachex2.append(x2)
cachey2.append(y2)
j +=1
teinschwing = cachet[0]
return plot(cachex2, cachey2, '.')

def poincare2(x1, y1):
if fabs(np.sin(x1)) <= 0.01 and y1 > 0:
return True
else:
return

Beim ausführen des Codes, haben wir für f1, doppelpendelz2 und für f2 doppelpendelz1 eingesetzt.
Es stekkt sich herraus, dass die Werte in die Höhe schießen und irgendwan so groß werden, und zu nan oder inf werden.

Wir sind euch sehr dankbar für Anregungen oder Erklärungen wie diese gekoppelten Bewegungsgleichungen miteinander zu handeln sind.
Vielen Dank im Vorraus!