Forum |  HardWare.fr | News | Articles | PC | Prix | S'identifier | S'inscrire | Aide | Shop Recherche
246 connectés 

  FORUM HardWare.fr
  Programmation
  Python

  Instabilité d'un système masse-ressort à 2 noeuds avec RK45

 



 Mot :   Pseudo :  
 
Bas de page
Auteur Sujet :

Instabilité d'un système masse-ressort à 2 noeuds avec RK45

n°2336598
Feitan21
Posté le 09-07-2019 à 15:23:07  profilanswer
 

:hello:  
 
Je cherche à modéliser des déformations d'un système avec un certain nombres de noeuds reliés entre eux par des ressorts.
 
J'ai un résultat un peu inattendu, lorsque mon noeud 0 est en position 0., 0., 0. et mon noeud 1 en position 1., 1., 1. tout fonctionne, et l'effet du ressort est bien cohérent. Le noeud 0 est toujours de masse infinie et le noeud 1 de masse 1 et j'applique une force similaire à la gravité en Y.
 
Le déplacement du noeud 1 en Y :  
https://reho.st/self/7f0473594ac08fbc57bf33261461dc6b3bd41b38.png
 
Par contre quand j'inverse la position des deux noeuds le ressort ne fait plus effet et mon noeud diverge vers l'infini.
 
Déplacement du noeud 1 en Y :
https://reho.st/self/2bc88ba7f871f77b0b1acdca62137da6e6df81bc.png
 
 
Je n'arrive pas à comprendre pourquoi j'obtiens ce résultat, je ne sais pas si c'est une instabilité numérique due au solver, ou si mon calcul de la force du ressort est erronée. J'ai pourtant essayé d'implementer une version standard avec des boucles imbriquées et une version vectorisée. Les deux donnent les mêmes résultats.
 
Voici mon code :  
 

Code :
  1. import numpy as np
  2. from scipy.spatial.distance import cdist
  3. import matplotlib.pyplot as plt
  4. import scipy.integrate
  5. X = np.array([[0., 0., 0.], # 0
  6.              [1., 1., 1.]]) # 1
  7. edges = np.array([[0, 1]])
  8. W = np.array([[0, 1],
  9.               [1, 0]])
  10. d = cdist(X, X, "euclidean" )
  11. mass = 1 * np.ones(X.shape[0])
  12. mass[0] = np.inf
  13. mass[1] = 1
  14. Force = np.zeros((X.shape[0], 3))
  15. Force[1, 1] = 9.81
  16. N = X.shape[0]
  17. dt = 0.1
  18. T = int(1e4)
  19. k = 5
  20. def compute_F_spring_vec(X, W, di):
  21.     d_inv = np.reciprocal(cdist(X, X, "euclidean" ))
  22.     d_inv[np.arange(np.shape(d_inv)[0]),np.arange(np.shape(d_inv)[0])] = 0
  23.     op1 = np.dot(W, X)
  24.     op2 = X * np.repeat(np.sum(W, axis=1), 3).reshape(X.shape)
  25.     op3 = np.dot(W * di * d_inv, X)
  26.     op4 = X * np.repeat(np.sum(W * di * d_inv, axis=1), 3).reshape(X.shape)
  27.     F_spring = - k * (-op1 + op2 + op3 - op4)
  28.     return F_spring
  29. def compute_F_spring(X, W, d):
  30.     F_spring = []
  31.     for i in range(X.shape[0]):
  32.         F_i = np.zeros(3)
  33.         for j in range(X.shape[0]):
  34.             if W[i,j]:
  35.                 F_i += -k * (np.linalg.norm(X[i, :] - X[j, :]) - d[i, j]) * (X[i, :] - X[j, :])/np.linalg.norm(X[i, :] - X[j, :])
  36.         F_spring.append(F_i)
  37.     return F_spring
  38. def fun(t, X):
  39.     X = X.reshape((int(X.shape[0]/3), 3))
  40.     F_spring = compute_F_spring_vec(X, W, d)
  41.     F = F_spring + Force
  42.     a = F / mass[:, np.newaxis]
  43.     X = X + a*t
  44.     return X.ravel()
  45. solver = scipy.integrate.RK45(fun, 0, X.ravel(), t_bound=100, rtol=0.001, vectorized=True)
  46. Ylast = np.inf
  47. X1 = []
  48. while solver.status == 'running' and solver.t < 10:
  49.     error = solver.step()
  50.     X1.append([solver.y[4], solver.t, solver.y[3], solver.y[5]])
  51.     if (np.abs(X1[-1][0] - Ylast))/np.abs(Ylast) < 1e-6:
  52.         break
  53.     else:
  54.         Ylast = X1[-1][0]


 
Quelques précisions :  
 
X contient les positions x,y,z de tous les nodes.  
Edges la liaison entre les nodes, ce qui me permet de déduire une matrice de connectivité W (1 les noeuds sont connectés, 0 ils ne le sont pas).
Et d la distance euclidienne entre chaque paire de noeuds
 
 
Pouvez-vous m'aider ?  
 
Merci


Message édité par Feitan21 le 09-07-2019 à 17:09:07
mood
Publicité
Posté le 09-07-2019 à 15:23:07  profilanswer
 

n°2336600
MaybeEijOr​Not
but someone at least
Posté le 09-07-2019 à 15:57:55  profilanswer
 

Sans lire le code, si tu places ta masse infinie en haut et que tu la soumets à une force de gravité, ça me semble logique que ça merde, puisse t-il y avoir un rapport ?


---------------
C'est en écrivant n'importe quoi qu'on devient n'importe qui.
n°2336609
Feitan21
Posté le 09-07-2019 à 17:07:06  profilanswer
 

MaybeEijOrNot a écrit :

Sans lire le code, si tu places ta masse infinie en haut et que tu la soumets à une force de gravité, ça me semble logique que ça merde, puisse t-il y avoir un rapport ?


 
Hello, merci pour ta réponse.
Je n'applique la force que sur le node 1, quelque-soit sa position. D'ailleurs si j'applique la gravité sur le noeud de masse infinie j'obtiens comme attendu un non déplacement puisque l'accélération est nulle.

n°2336619
MaybeEijOr​Not
but someone at least
Posté le 09-07-2019 à 19:56:54  profilanswer
 

Pour inverser les positions tu fais quoi ?

 

X = np.array([[1., 1., 1.], # 0
             [0., 0., 0.]]) # 1


Et :

Code :
  1. mass[0] = 1
  2. mass[1] = np.inf
 

?

 

EDIT : et qu'est-ce que ça fait si tu fais :

Code :
  1. X = np.array([[0., 0., 0.], # 0
  2.              [-1., -1., -1.]]) # 1

Message cité 1 fois
Message édité par MaybeEijOrNot le 09-07-2019 à 20:04:32

---------------
C'est en écrivant n'importe quoi qu'on devient n'importe qui.
n°2336658
Feitan21
Posté le 10-07-2019 à 14:41:48  profilanswer
 

MaybeEijOrNot a écrit :

Pour inverser les positions tu fais quoi ?
 

X = np.array([[1., 1., 1.], # 0
             [0., 0., 0.]]) # 1


Et :

Code :
  1. mass[0] = 1
  2. mass[1] = np.inf


 
?
 
EDIT : et qu'est-ce que ça fait si tu fais :

Code :
  1. X = np.array([[0., 0., 0.], # 0
  2.              [-1., -1., -1.]]) # 1



 
Hello,  
 
Pour inverser les positions je change juste les coordonnées x,y,z de mon noeud 0 et 1  
 

X = np.array([[1., 1., 1.], # 0
             [0., 0., 0.]]) # 1


 
Je ne touche pas à la masse.
 
D'ailleurs si je rentre tes paramètres de masse, j'ai bien mon point 1 fixe, et mon point 0 qui est cohérent.
 
 
Si je teste  
 

Code :
  1. X = np.array([[0., 0., 0.], # 0
  2.              [-1., -1., -1.]]) # 1


 
J'ai bien convergence, et cette valeur est bien cohérente puisque ma force est positive. Je tombe sur le même résultat qu'avec un noeud 1 placé en 1, 1, 1
 
https://reho.st/thumb/self/4e61eb0085ab566bfddbc69506850ad20130479d.png


Message édité par Feitan21 le 10-07-2019 à 14:43:41
n°2336663
MaybeEijOr​Not
but someone at least
Posté le 10-07-2019 à 15:49:19  profilanswer
 

Et tu obtiens quoi comme matrice d'accélération dans les deux configurations ?


---------------
C'est en écrivant n'importe quoi qu'on devient n'importe qui.
n°2336669
Feitan21
Posté le 10-07-2019 à 16:24:47  profilanswer
 

Je sors les images de la matrice d'accélération en Y pour le noeud 1 uniquement, le noeud 0 étant fixe dans les deux cas:  
 
Pour la bonne partie, celle qui converge :  
https://reho.st/self/6f21c8b3aff861169ef341499211e28e385f4e95.png
 
Et pour celle qui diverge :  
 
https://reho.st/self/4acfbaae146eeb85bbb661e582a6d31e66ce14ad.png
 
En zoomant on remarque que le début est correct mais l'instabilité augmente au fur et à mesure avec des oscillations :  
 
https://reho.st/self/80e6de0cec583df62e697f1471a9a941d7f08b4b.png
 
 
Et les données brutes, avec la matrice d'accélération en x,y,z pour le noeud 1 dans les deux cas : https://docs.google.com/spreadsheet [...] sp=sharing

n°2336673
MaybeEijOr​Not
but someone at least
Posté le 10-07-2019 à 18:41:01  profilanswer
 

Si ce n'est pas la masse c'est peut-être la force.
 
Si tu inverses ta configuration, ta force ne devrait pas être :

Code :
  1. Force[0, 1] = 9.81


 :??:


---------------
C'est en écrivant n'importe quoi qu'on devient n'importe qui.
n°2336676
Feitan21
Posté le 10-07-2019 à 19:06:33  profilanswer
 

MaybeEijOrNot a écrit :

Si ce n'est pas la masse c'est peut-être la force.
 
Si tu inverses ta configuration, ta force ne devrait pas être :

Code :
  1. Force[0, 1] = 9.81


 :??:


 
Non ma force est exacte, je ne change pas la position des noeuds dans ma liste, juste mes positions. C'est ce qui me permet d'appliquer toujours ma force en [1,1], soit la deuxième dimension du node 1. le node 0 étant toujours fixe. De plus j'ai bien mon noeud 0 qui est fixe, donc je ne pense pas que ce soit un problème d'application de la force.
 
Après pour une erreur dans le calcul de la force c'est possible, mais je trouve ça étonnant que ça fonctionne dans de nombreuses conditions et qu'il y ait quelques configurations ou ça ne fonctionne pas.


Message édité par Feitan21 le 10-07-2019 à 19:09:43

Aller à :
Ajouter une réponse
  FORUM HardWare.fr
  Programmation
  Python

  Instabilité d'un système masse-ressort à 2 noeuds avec RK45

 

Sujets relatifs
XSLT 1.0 - Grouper par liste de noeuds identiquesSystème d'achievements/badges et statistiques
Système de log de champs de BD mis à jourphp - script de récupération infos système
Programme simple pour récupérer informations système d'une machine linSystème de recommandation sur hébergeur mutualisé
Envoi d'emails en masseProjet de système à développer
Systeme de sessionsystème référentiel / masse vitesse et position relative.
Plus de sujets relatifs à : Instabilité d'un système masse-ressort à 2 noeuds avec RK45


Copyright © 1997-2018 Hardware.fr SARL (Signaler un contenu illicite) / Groupe LDLC / Shop HFR