Système couplé de 4 équations différentielles - Python

Question

J'ai le système couplé de 4 équations différentielles dans l'image. J'ai 4 fonctions (xG; yG; gamma; beta) et leurs dérivés. Ils sont tous fonction de la même variable indépendante t.

J'essaye de le résoudre avec odeint. Le problème est que, pour ce faire, je pense que je dois exprimer le système de manière à ce que chaque dérivée seconde ne dépende pas d'autres dérivées secondes. Cela implique une quantité de mathématiques qui me mènera sûrement à une erreur quelque part (j'ai essayé!).

Savez-vous comment je pourrais soit:

1. Résoudre ce système d'équations différentielles tel quel?
2. ou demander à python d'isoler les dérivés secondaires pour moi?

Je joins mon code de test

Merci

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

Vous devez réécrire toutes vos dérivées du second ordre en tant que dérivées du premier ordre et résoudre 8 ODE ensemble:

Ensuite, vous avez besoin des conditions initiales pour tous les dérivés, mais il semble que vous les ayez déjà. Pour info, votre code ne fonctionne pas ( line 71: NameError: name 'Xg2' is not defined), veuillez le vérifier.

Pour plus d'informations, reportez-vous également à la résolution numérique de l'ODE d'ordre 2 .

EDIT # 1: À la première étape, vous devez découpler le système d'équations. Bien que vous puissiez le résoudre manuellement, je ne le recommanderais pas, alors utilisons le sympymodule:

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

maintenant nous avons tous les dérivés du 2ème ordre découplés. Par exemple, l'expression pour beta2est

donc il (et tous les autres dérivés du 2ème ordre aussi) a la forme

notez qu'il n'y a pas de dépendance sur xgou yg.

Introduisons les deux nouvelles variables, bet k:

alors le

devient

et le système complet d'EDO à résoudre est

Désormais, tous les ODE dépendent de quatre variables qui ne sont dérivées de rien. Aussi, puisque les xget ygsont dégénérés, il n'y a aussi que 6 équations au lieu de 8. Cependant, on peut réécrire ces deux équations de la même manière que gammaet betapour obtenir le système complet de 8 équations, et l'intégrer ensemble.