Python : système d'équations différentielles liées

[cmole]

Bonjour à tous,

Je suis (une nouvelle fois) en train d’essayer de résoudre un système d’équations différentielles à l’aide de Python, et je dois avouer que je sèche un peu… J’ai utilisé ChatGPT, et même si mes compétences en la matière sont très limitées, je ne suis pas convaincu de la piste qui m’a été donnée (piste qui diffère si je repose la question dans différentes conversations, qui plus est).

Je vais essayer de vous décrire mon problème du mieux que je peux, et vous donner le code que j’ai utilisé jusque-là. Voici le système d’équations que j’ai à résoudre :
N1 = A.x’’cos(x) + B.x’²sin(x) + C
N2 = N1.D.sin(x) + E.x’²
x’’ = F.N2 + G.N1.cos(x) + H.x’² + I

Avec x(t) la fonction inconnue (on pourra dire que x est une position, x’ une vitesse et x’’ une accélération) ; A, B, … , I sont des constantes, et N1 et N2 sont des paramètres variables, dépendants de x(t) et ses dérivées.
Je m’y prends surement mal, mais je n’arrive pas à résoudre ce problème, et je n’arrive pas à trouver de « documentation » me satisfaisant pour me guider. Pour me simplifier la tâche j’ai essayé de décomposer mon problème : j’ai donc supposé que F = 0, comme ça je peux ignorer N2, puis « à la main » j’ai injecté N1 dans l’expression de x’’, puis j’ai résolu l’équation à l’aide de Python, aucun problème jusque-là. Ensuite j’ai voulu résoudre ce système de deux équations (N1 et x’’) pour voir si je trouvais le même résultat que lors de ma première résolution, et là ça cloche. Avec la résolution de ChatGPT j’ai une erreur supérieure à 5%, alors certes ce n’est pas énorme, mais je n’ai pas encore intégré N2 ! Voici le code :

# Définition du système d'équation différentielle
def systeme(t, X):
x1, x2 = X

def compute_N1 (x1, x2, dx2_dt):
return(A*dx2_dt*cos(x1) + B*x2²*sin(x) + C

dx2_dt_guess = 0.0

for _ in range(10):
N1 = compute_N1(x1, x2, dx2_dt_guess)
dx2_dt_guess = G*N1*cos(x1) + H*x2² + I

dx1_dt = x2
dx2_dt = dx2_dt_guess
return [dx1_dt, dx2_dt]

Ensuite résolution avec solve_ivp, mais ce n’est pas très important.
Ma question est donc : la méthode employée vous semble-t-elle cohérente ? si non, que feriez-vous ?

Merci d’avance !

NB: Visiblement la tabulation ne se fait pas, j'espère que le code reste clair, désolé...
-----

[gts2]

Bonjour,

Vous avez un problème un peu tordu...

Deux solutions
- vous résolvez à la main par une méthode d'Euler ou plus évolué et comme vous contrôlez le détail des calculs, vous n'aurez aucun mal à introduire N1 et N2
- vous utilisez scipy ou numpy, dans ce cas une possibilité est de considérer
N1 = A x2' \cos(x1) + B x2^2\sin(x1) + C
N2 = N_1 D \sin(x1) + E x2^2
x2'= F N2 + G N1 \cos(x) + H x2^2 + I
comme un système linéaire en (N1,N2,x2'), résoudre ce système et renvoyer x2,x2'

[polo974]

Pour la lisibilité il faut passer en mode édition avancée et utiliser la balise cd code (bouton #)
ce qui donne:

Bonjour à tous,
snip snip snip ...
Voici le code :

# Définition du système d'équation différentielle
def systeme(t, X):
            x1, x2 = X

        def compute_N1 (x1, x2, dx2_dt):
                return(A*dx2_dt*cos(x1) + B*x2²*sin(x) + C

      dx2_dt_guess = 0.0

       for _ in range(10):
              N1 = compute_N1(x1, x2, dx2_dt_guess)
                dx2_dt_guess = G*N1*cos(x1) + H*x2² + I
    
    dx1_dt = x2
    dx2_dt = dx2_dt_guess
    return [dx1_dt, dx2_dt]

Ensuite résolution avec solve_ivp, mais ce n’est pas très important.
Ma question est donc : la méthode employée vous semble-t-elle cohérente ? si non, que feriez-vous ?

Merci d’avance !

NB: Visiblement la tabulation ne se fait pas, j'espère que le code reste clair, désolé...

Edit: Mais là, l'indentation est un peu pourrie...

[cmole]

Pour la lisibilité il faut passer en mode édition avancée et utiliser la balise cd code (bouton #)
ce qui donne:

Je le saurai pour la prochaine fois, merci !

- vous utilisez scipy ou numpy, dans ce cas une possibilité est de considérer
N1 = A x2' \cos(x1) + B x2^2\sin(x1) + C
N2 = N_1 D \sin(x1) + E x2^2
x2'= F N2 + G N1 \cos(x) + H x2^2 + I
comme un système linéaire en (N1,N2,x2'), résoudre ce système et renvoyer x2,x2'

C'est exactement cette méthode que je souhaiterais suivre, auriez-vous un exemple ou une piste pour m'aider à implémenter ça ?

[A voir en vidéo sur Futura]

[MissJenny]

quand j'ai eu à faire un travail similaire j'ai utilisé la bibliiothèque "Numerical Recipes". C'est écrit en c, je ne sais pas comment on fait pour appeler ces fonctions depuis Python mais ça doit être possible.

[gts2]

Bonjour,

Etape 1 (mathématique) écrire les équations sous la forme AX=B qqch du genre :

; ;

Calcul bien sûr à vérifier ...

Etape 2 (Python) utilisez numpy.linalg.solve

\begin{lstlisting}
a = np.array([[ , , ], [ , , ],[ , , ]])
b = np.array([ , , ])
x = np.linalg.solve(a, b)
dx2_dt = x[2]
\end{lstlisting}

[Biname]

Salut,
L'IA(deepseek ici) donne une solution numérique sous python, plot et régression en 30 secondes avec un copier-coller de l'énoncé.

 Cliquez pour afficher

import numpy as np
from scipy.integrate import solve_ivp
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# Constantes
A, B, C, D, E, F, G, H, I = 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  # Remplace par les vraies valeurs

# Définir le système d'équations
def system(t, y):
    y1, y2 = y  # y1 = x, y2 = x'
    
    # Calculer N1
    N1 = A * y2 * np.cos(y1) + B * y2**2 * np.sin(y1) + C
    
    # Calculer N2
    N2 = N1 * D * np.sin(y1) + E * y2**2
    
    # Calculer x'' (y2')
    y2_prime = F * N2 + G * N1 * np.cos(y1) + H * y2**2 + I
    
    # Retourner les dérivées
    return [y2, y2_prime]

# Conditions initiales
y0 = [0.0, 1.0]  # x(0) = 0, x'(0) = 1

# Intervalle de temps
t_span = (0, 10)  # De t=0 à t=10

# Résolution
sol = solve_ivp(system, t_span, y0, method='RK45', t_eval=np.linspace(0, 10, 1000))

# Ajuster un polynôme à x(t)
def polynomial(t, *coeffs):
    return np.polyval(coeffs, t)

# Degré du polynôme (ajuste selon tes besoins)
degree = 5
p0 = np.ones(degree + 1)  # Valeurs initiales pour les coefficients

# Ajustement
popt, _ = curve_fit(polynomial, sol.t, sol.y[0], p0=p0)

# Afficher les coefficients du polynôme
print("Coefficients du polynôme ajusté :", popt)

# Tracer la solution numérique et l'approximation polynomiale
plt.plot(sol.t, sol.y[0], label='x(t) (numérique)')
plt.plot(sol.t, polynomial(sol.t, *popt), '--', label='x(t) (polynomial)')
plt.xlabel('Temps')
plt.ylabel('x(t)')
plt.legend()
plt.show()
print (F"\n\n     *** FINI ***\n\n")

Code solution exacte (trop long ou ne trouve pas) :

 Cliquez pour afficher

import sympy as sp

# Définir les symboles
t = sp.symbols('t')
x = sp.Function('x')(t)
A, B, C, D, E, F, G, H, I = sp.symbols('A B C D E F G H I')

# Définir les dérivées
x_dot = sp.diff(x, t)
x_ddot = sp.diff(x_dot, t)

# Définir N1 et N2
N1 = A * x_ddot * sp.cos(x) + B * x_dot**2 * sp.sin(x) + C
N2 = N1 * D * sp.sin(x) + E * x_dot**2

# Définir l'équation différentielle pour x''
eq = sp.Eq(x_ddot, F * N2 + G * N1 * sp.cos(x) + H * x_dot**2 + I)

# Afficher l'équation
print("Équation différentielle :")
sp.pprint(eq)

# Essayer de résoudre symboliquement (peut échouer)
try:
    solution = sp.dsolve(eq, x)
    print("Solution symbolique :")
    sp.pprint(solution)
except Exception as e:
    print("Une solution symbolique n'a pas pu être trouvée :", e)

[gts2]

Le problème de la solution proposée est qu'elle ne correspond pas à la question posée :

question : N1 = A.x’’.cos(x) + B.x’²sin(x) + C
solution : N1 = A.x’.cos(x) + B.x’²sin(x) + C

[Biname]

Oui, il a triché . Amusant : Mistral et 4o font la même erreur.

[pm42]

Le problème de la solution proposée est qu'elle ne correspond pas à la question posée :

Normal : poser ce genre de question à un LLM qui n'est pas fait pour ça au lieu d'un LRM qui a de meilleures chances de faire quelque chose de correct.
Ce qui est curieux, c'est que cela a été expliqué de nombreuses fois mais on récupère quand même ces copier-coller sans valeur ajoutée.

[gts2]

Bonjour,

quand j'ai eu à faire un travail similaire j'ai utilisé la bibliiothèque "Numerical Recipes". C'est écrit en c, je ne sais pas comment on fait pour appeler ces fonctions depuis Python mais ça doit être possible.

En fait le titre est trompeur : il ne s'agit pas de "système d'équations différentielles" mais d'UNE équation différentielle donnée par un système.

Sinon pour Python, c'est la même procédure pour les systèmes d'équation ou une équation : dans un cas on passe comme argument des scalaires, dans l'autre des vecteurs : Python n'est pas un langage typé.

[Biname]

Salut,
L'IA a raison, on a juste pas compris le changement de variables

[gts2]

l'IA a peut-être changé d'avis entre temps mais la solution initialement proposée est :

def system(t, y):
y1, y2 = y # y1 = x, y2 = x'
# Calculer N1
N1 = A * y2 * np.cos(yl) + B * y2**2 * np.sin(yl) + C

Donc clairement différent de ce qui est proposé ci-dessus.

[cmole]

Bonjour à tous,

Un grand merci pour tous pour vos retours et votre aide. Je vais me repencher un peu sur mon problème ce weekend, et ne manquerai pas de vous partager mon avancée !

[Biname]

l'IA a peut-être changé d'avis entre temps mais la solution initialement proposée est :

def system(t, y):
y1, y2 = y # y1 = x, y2 = x'
# Calculer N1
N1 = A * y2 * np.cos(yl) + B * y2**2 * np.sin(yl) + C

Donc clairement différent de ce qui est proposé ci-dessus.

Oui. S'il me reste un chapeau, je vais le manger !

[Biname]

Salut,
Cette fois-ci, je lui ai demandé de réfléchir .
Il a sorti x'' = y2' = F(A,..,I, y1) qui est une fraction dont D peut valoir 0. Je lui a demandé de faire l'étude pour D != 0
On voit immédiatement que pour A = 0, D > 0, ce qui annule x'' dans N1=

Soyons prudent : sauf erreursssss :

 Cliquez pour afficher

Enoncé



Variables



Réécriture



Condition et raisonnement




Le Code Python

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Define constants (replace with actual values)
A, B, C, D, E, F, G, H, I = 1, 1, 1, 1, 1, 1, 1, 1, 1

# Define the system of first-order ODEs
def system(t, y):
    y1, y2 = y
    numerator = (F * B * D * y2**2 * np.sin(y1)**2 +
                 F * C * D * np.sin(y1) +
                 F * E * y2**2 +
                 G * B * y2**2 * np.sin(y1) * np.cos(y1) +
                 G * C * np.cos(y1) +
                 H * y2**2 + I)
    denominator = 1 - F * A * D * np.cos(y1) * np.sin(y1) - G * A * np.cos(y1)**2
    dy1dt = y2
    dy2dt = numerator / denominator
    return [dy1dt, dy2dt]

# Set initial conditions
y0 = [0, 0]  # x(0) = 0, x'(0) = 0

# Set time span for integration
t_span = (0, 10)

# Solve the ODE system
solution = solve_ivp(system, t_span, y0, dense_output=True)

# Plot the results
t = np.linspace(t_span[0], t_span[1], 300)
y = solution.sol(t)
plt.figure(figsize=(10, 6))
plt.plot(t, y[0], label='x(t)')
plt.plot(t, y[1], label="x'(t)")
plt.xlabel('Time t')
plt.ylabel('Values')
plt.title('Solution of the ODE System')
plt.legend()
plt.grid(True)
plt.show()

[JPL]

Il serait peut-être temps de laisser des neurone humains réfléchir ! Stop à l’IA maintenant

[cmole]

Bonsoir,

Je voulais tous vous remercier pour votre aide sur ce sujet, j'ai réussi à résoudre mon problème ! J'ai suivi la méthode indiquée par gts2, et je dois dire que sa remarque

il ne s'agit pas de "système d'équations différentielles" mais d'UNE équation différentielle donnée par un système.

est plus que pertinente. Je m'en veux un peu de ne pas y avoir pensé tant ça me semble "logique" maintenant, j'y ferais attention la prochaine fois.

Encore merci !