Sprzężony układ 4 równań różniczkowych - Python

Question

Na obrazku mam sprzężony układ 4 równań różniczkowych. Mam 4 funkcje (xG; yG; gamma; beta) i ich pochodne. Wszystkie są funkcją tej samej zmiennej niezależnej t.

Próbuję to rozwiązać za pomocą odeint. Problem polega na tym, że aby to zrobić, myślę, że muszę wyrazić system w taki sposób, aby każda druga pochodna nie była zależna od innych drugich pochodnych. Wiąże się to z dużą ilością matematyki, która z pewnością doprowadzi mnie do błędu (próbowałem!).

Czy wiesz, jak mogłem:

1. Jak rozwiązać ten układ równań różniczkowych?
2. lub poprosić Pythona o wyodrębnienie dla mnie drugiej pochodnej?

Załączam kod testowy

Dzięki

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

Musisz przepisać wszystkie pochodne drugiego rzędu na pochodne pierwszego rzędu i rozwiązać razem 8 ODE:

Następnie potrzebujesz warunków początkowych dla wszystkich instrumentów pochodnych, ale wydaje się, że już je masz. FYI, Twój kod nie działa ( line 71: NameError: name 'Xg2' is not defined), sprawdź to.

Aby uzyskać więcej informacji, zobacz numeryczne rozwiązywanie ODE drugiego rzędu .

EDYCJA # 1: W pierwszym kroku musisz oddzielić układ równań. Chociaż mógłbyś rozwiązać to ręcznie, nie polecam, więc użyjmy sympymodułu:

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])

teraz mamy odsprzęgnięte wszystkie pochodne drugiego rzędu. Na przykład wyrażenie for beta2is

więc to (i wszystkie inne pochodne drugiego rzędu również) ma postać

zauważ, że nie ma zależności od xglub yg.

Przedstawmy dwie nowe zmienne bi k:

wtedy

staje się

a cały system ODE do rozwiązania jest

Teraz wszystkie ODE są zależne od czterech zmiennych, które nie są pochodnymi niczego. Ponadto, ponieważ xgi ygsą zdegenerowane, istnieje również tylko 6 równań zamiast 8. Jednak można przepisać te dwa równania w taki sam sposób, jak gammai betauzyskać pełny układ 8 równań i scalić je razem.