Gekoppeltes System von 4 Differentialgleichungen - Python

Question

Ich habe das gekoppelte System von 4 Differentialgleichungen im Bild. Ich habe 4 Funktionen (xG; yG; Gamma; Beta) und deren Derivate. Sie sind alle Funktionen derselben unabhängigen Variablen t.

Ich versuche es mit Odeint zu lösen. Das Problem ist, dass ich dazu das System so ausdrücken muss, dass jede zweite Ableitung nicht von anderen zweiten Ableitungen abhängig ist. Dies beinhaltet eine Menge Mathe, die mich sicher irgendwo zu einem Fehler führen wird (ich habe es versucht!).

Weißt du, wie ich entweder könnte:

1. Dieses Differentialgleichungssystem so lösen, wie es ist?
2. oder Python dazu bringen, die zweiten Ableitungen für mich zu isolieren?

Ich füge meinen Testcode hinzu

Vielen Dank

import numpy
import math
from numpy import loadtxt
from pylab import figure,  savefig
import matplotlib.pyplot as plt
# Use ODEINT to solve the differential equations defined by the vector field
from scipy.integrate import odeint



def vectorfield(w, t, p):
    """
    Defines the differential equations for the coupled system.

    Arguments:
        w :  vector of the state variables:
                  w = [Xg, Xg1 Yg, Yg1, Gamma, Gamma1, Beta, Beta1]
        t :  time
        p :  vector of the parameters:
                  p = [m, rAG, Ig,lcavo]
    """
#Xg is position ; Xg1 is the first derivative ; Xg2 is the second derivative (the same for the other functions)
        Xg, Xg1,  Yg, Yg1, Gamma, Gamma1, Beta, Beta1 = w
        Xg2=-(Ig*Gamma2*math.cos(Beta))/(rAG*m*(-math.cos(Gamma)*math.sin(Beta)+math.sin(Gamma)*math.cos(Beta)))
        Yg2=-(Ig*Gamma2*math.sin(Beta))/(rAG*m*(-math.cos(Gamma)*math.sin(Beta)+math.sin(Gamma)*math.cos(Beta)))-9.81
        Gamma2=((Beta2*lcavo*math.sin(Beta))+(Beta1**2*lcavo*math.cos(Beta))+(Xg2)-(Gamma1**2*rAG*math.cos(Gamma)))/(rAG*math.sin(Gamma))
        Beta2=((Yg2)+(Gamma2*rAG*math.cos(Gamma))-(Gamma1**2*rAG*math.sin(Gamma))+(Beta1**2*lcavo*math.sin(Beta)))/(lcavo*math.cos(Beta))
        m, rAG, Ig,lcavo, Xg2,  Yg2, Gamma2, Beta2 = p
    
    
    # Create f = (Xg', Xg1' Yg', Yg1', Gamma', Gamma1', Beta', Beta1'):
    f = [Xg1,
         Xg2,
         Yg1, 
         Yg2, 
         Gamma1, 
         Gamma2, 
         Beta1, 
         Beta2]
         
    return f

    


# Parameter values
m=2.722*10**4
rAG=2.622
Ig=3.582*10**5
lcavo=4
# Initial conditions
Xg = 0.0
Xg1 = 0
Yg = 0.0
Yg1 = 0.0
Gamma=-2.52
Gamma1=0
Beta=4.7
Beta1=0

# ODE solver parameters
abserr = 1.0e-8
relerr = 1.0e-6
stoptime = 5.0
numpoints = 250

#create the time values
t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]
Deltat=t[1]
# Pack up the parameters and initial conditions:
p = [m, rAG, Ig,lcavo, Xg2,  Yg2, Gamma2, Beta2]
w0 = [Xg, Xg1,  Yg, Yg1, Gamma, Gamma1, Beta, Beta1]

# Call the ODE solver.
wsol = odeint(vectorfield, w0, t, args=(p,),
              atol=abserr, rtol=relerr)

Answer 1

Sie müssen alle Ihre Derivate zweiter Ordnung als solche erster Ordnung umschreiben und 8 ODE zusammen lösen:

Dann brauchen Sie Anfangsbedingungen für alle Derivate, aber es scheint, dass Sie bereits haben. Zu Ihrer Information, Ihr Code läuft nicht ( line 71: NameError: name 'Xg2' is not defined), bitte überprüfen Sie es.

Weitere Informationen finden Sie unter numerisches Lösen der ODE 2. Ordnung .

EDIT # 1: Im ersten Schritt müssen Sie das Gleichungssystem entkoppeln. Sie könnten es zwar manuell lösen, aber ich würde es nicht empfehlen. Verwenden wir also das sympyModul:

import sympy as sm
from sympy import symbols

# define symbols. I assume all the variables are real-valued, this helps the solver. If not, I believe the result will be the same, but just calculated slower
Ig, gamma, gamma1, gamma2, r, m, beta, beta1, beta2, xg2, yg2, g, l = symbols('I_g, gamma, gamma1, gamma2, r, m, beta, beta1, beta2, xg2, yg2, g, l', real = True)

# define left hand sides as expressions
# 2nd deriv of gamma
g2 = (beta2 * l * sm.sin(beta) + beta1**2 *l *sm.cos(beta) + xg2 - gamma1**2 *r * sm.cos(gamma))/(r*sm.sin(gamma))
# 2nd deriv of beta
b2 = (yg2 + gamma2 * r * sm.cos(gamma) - gamma1**2 *r * sm.sin(gamma) + beta1**2 *l *sm.sin(beta))/(l*sm.cos(beta))
# 2nd deriv of xg
x2 = -Ig*gamma2*sm.cos(beta)/(r*m*(-sm.sin(beta)*sm.cos(gamma) + sm.sin(gamma)*sm.cos(beta)))
# 2nd deriv of yg
y2 = -Ig*gamma2*sm.sin(beta)/(r*m*(-sm.sin(beta)*sm.cos(gamma) + sm.sin(gamma)*sm.cos(beta))) - g

# now let's solve the system of four equations to decouple second order derivs
# gamma2 - g2 means "gamma2 - g2 = 0" to the solver. The g2 contains gamma2 by definition
# one could define these equations the other way, but I prefer this form
result = sm.solve([gamma2-g2,beta2-b2,xg2-x2,yg2-y2],
                  # this line tells the solver what variables we want to solve to
                  [gamma2,beta2,xg2,yg2] )
# print the result
# note that it is long and ugly, but you can copy-paste it as python code
for res in result:
    print(res, result[res])

Jetzt haben wir alle Derivate 2. Ordnung entkoppelt. Zum Beispiel für den Ausdruck beta2ist

so hat es (und alle anderen Derivate 2. Ordnung auch) die Form

Beachten Sie, dass keine Abhängigkeit von xgoder besteht yg.

Lassen Sie uns die beiden neuen Variablen vorstellen bund k:

dann

wird das

und das vollständige System der zu lösenden ODEs ist

Jetzt sind alle ODEs von vier Variablen abhängig, die keine Ableitungen von irgendetwas sind. Da die xgund ygentartet sind, gibt es auch nur 6 Gleichungen anstelle von 8. Man kann diese beiden Gleichungen jedoch auf die gleiche Weise wie gammaund umschreiben beta, um ein vollständiges System von 8 Gleichungen zu erhalten und sie zusammen zu integrieren.