Ressorts
Description
voici un petit programme qui montre graphiquement la solution du probleme suivant :
le probleme est de resoudre l'equation differentielle d'un systeme compose de 3 masses A, B, et C dont chacune est reliees par un ressort (ki,l0i) avec i=A,B ou C, et k la constante de raideur et l0 la longeur a vide
je mets ci-dessous le fonction de l'equation differentielle (il faut voir la suite du code dans main.c pour les conditions initiales)
Les equa. diff. sont resolues avec du Runge-Kutta a pas adaptatif.
le probleme est de resoudre l'equation differentielle d'un systeme compose de 3 masses A, B, et C dont chacune est reliees par un ressort (ki,l0i) avec i=A,B ou C, et k la constante de raideur et l0 la longeur a vide
je mets ci-dessous le fonction de l'equation differentielle (il faut voir la suite du code dans main.c pour les conditions initiales)
Les equa. diff. sont resolues avec du Runge-Kutta a pas adaptatif.
Source / Exemple :
//-------------------------------------------------
typedef struct tagRESSORT
{
double k; // constante de raideur
double l0; // longeur a vide
}RESSORT,*P_RESSORT,**PP_RESSORT;
// il y a 3 masses A, B, et C
// liees par des ressorts
typedef struct tagEQUATION
{
// les donnees
double mA,mB,mC; // masses
RESSORT rAB,rAC,rBC; // les trois ressorts
// informations pour le dessin
COLORREF cA,cB,cC; // couleurs des masses
COLORREF cRessorts; // couleurs des ressorts
DWORD t0;
double alpha; // facteur d'amortissement
// l'integrateur d'equations differentielles
ODE_RK rk;
}EQUATION,*P_EQUATION,**PP_EQUATION;
//-------------------------------------------------
// 12 valeurs : (3 masses) * (position (x,y) + vitesse (vx,vy)) = (3)*(2+2) = 12
//-------------------------------------------------
#define INDEX_XA 0
#define INDEX_YA 1
#define INDEX_VXA 2
#define INDEX_VYA 3
#define INDEX_XB 4
#define INDEX_YB 5
#define INDEX_VXB 6
#define INDEX_VYB 7
#define INDEX_XC 8
#define INDEX_YC 9
#define INDEX_VXC 10
#define INDEX_VYC 11
//-------------------------------------------------
// 12 valeurs : (3 masses) * (vitesse (vx,vy) + acceleration (ax,ay)) = (3)*(2+2) = 12
//-------------------------------------------------
#define DINDEX_VXA 0
#define DINDEX_VYA 1
#define DINDEX_AXA 2
#define DINDEX_AYA 3
#define DINDEX_VXB 4
#define DINDEX_VYB 5
#define DINDEX_AXB 6
#define DINDEX_AYB 7
#define DINDEX_VXC 8
#define DINDEX_VYC 9
#define DINDEX_AXC 10
#define DINDEX_AYC 11
//-------------------------------------------------
// on cherche a resoudre Y'=f(X,Y)
// avec X un scalaire, et Y et Y' sont des vecteurs (12 dimensions) :
// Y = ( xA ) Y' = ( vxA ) <-- masse A
// ( yA ) ( vyA )
// ( vxA ) ( axA )
// ( vyA ) ( ayA )
// ( xB ) ( vxB ) <-- masse B
// ( yB ) ( vyB )
// ( vxB ) ( axB )
// ( vyB ) ( ayB )
// ( xC ) ( vxC ) <-- masse C
// ( yC ) ( vyC )
// ( vxC ) ( axC )
// ( vyC ) ( ayC )
//
// ici dans cette fonction <res> est le vecteur Y' que l'on doit calculer
BOOL funcEquaDiff(P_VECTND res,double x,P_VECTND y,void *param)
{
P_EQUATION eq;
double FAB,FAC,FBC; // les forces (en norme relative +/-)
double lAB,lAC,lBC; // les longueurs
double angleAB,angleAC,angleBC; // les angles des droites (AB), (AC), et (BC)
double dxAB,dyAB; // variables utiles (calcul des longueurs et des angles)
double dxAC,dyAC; // variables utiles (calcul des longueurs et des angles)
double dxBC,dyBC; // variables utiles (calcul des longueurs et des angles)
eq = (P_EQUATION)param;
AssertPointer(eq);
// le probleme des vitesses est trivial,
// puisque l'on nous les donne
res->coord[DINDEX_VXA] = y->coord[INDEX_VXA]; // masse A
res->coord[DINDEX_VYA] = y->coord[INDEX_VYA];
res->coord[DINDEX_VXB] = y->coord[INDEX_VXB]; // masse B
res->coord[DINDEX_VYB] = y->coord[INDEX_VYB];
res->coord[DINDEX_VXC] = y->coord[INDEX_VXC]; // masse C
res->coord[DINDEX_VYC] = y->coord[INDEX_VYC];
// passons aux choses serieuses, l'accelerations des masses
// est donnees par le PFD (principe fondamental de la dynamique)
// a=F/m, de plus on sait calculer F pour des ressorts : F = - k.(l-l0)
// * on calcule les longueurs
dxAB = y->coord[INDEX_XB] - y->coord[INDEX_XA]; // longueur entre A et B
dyAB = y->coord[INDEX_YB] - y->coord[INDEX_YA];
lAB = sqrt(dxAB*dxAB + dyAB*dyAB);
dxAC = y->coord[INDEX_XC] - y->coord[INDEX_XA]; // longueur entre A et C
dyAC = y->coord[INDEX_YC] - y->coord[INDEX_YA];
lAC = sqrt(dxAC*dxAC + dyAC*dyAC);
dxBC = y->coord[INDEX_XC] - y->coord[INDEX_XB]; // longueur entre B et C
dyBC = y->coord[INDEX_YC] - y->coord[INDEX_YB];
lBC = sqrt(dxBC*dxBC + dyBC*dyBC);
// * on calcule les forces
FAB = - eq->rAB.k * (lAB - eq->rAB.l0);
FAC = - eq->rAC.k * (lAC - eq->rAC.l0);
FBC = - eq->rBC.k * (lBC - eq->rBC.l0);
// * on calcule les angles (sens trigonometrique)
angleAB = atan2(dyAB,dxAB);
angleAC = atan2(dyAC,dxAC);
angleBC = atan2(dyBC,dxBC);
// * on calcule les accelerations (attention aux signes)
res->coord[DINDEX_AXA] = ((- FAB*cos(angleAB) - FAC*cos(angleAC))/eq->mA) - eq->alpha*y->coord[INDEX_VXA]; // masse A
res->coord[DINDEX_AYA] = ((- FAB*sin(angleAB) - FAC*sin(angleAC))/eq->mA) - eq->alpha*y->coord[INDEX_VYA];
res->coord[DINDEX_AXB] = ((+ FAB*cos(angleAB) - FBC*cos(angleBC))/eq->mB) - eq->alpha*y->coord[INDEX_VXB]; // masse B
res->coord[DINDEX_AYB] = ((+ FAB*sin(angleAB) - FBC*sin(angleBC))/eq->mB) - eq->alpha*y->coord[INDEX_VYB];
res->coord[DINDEX_AXC] = ((+ FAC*cos(angleAC) + FBC*cos(angleBC))/eq->mC) - eq->alpha*y->coord[INDEX_VXC]; // masse C
res->coord[DINDEX_AYC] = ((+ FAC*sin(angleAC) + FBC*sin(angleBC))/eq->mC) - eq->alpha*y->coord[INDEX_VYC];
return TRUE;
} // funcEquaDiff()
Vous n'êtes pas encore membre ?
inscrivez-vous, c'est gratuit et ça prend moins d'une minute !
Les membres obtiennent plus de réponses que les utilisateurs anonymes.
Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.
Le fait d'être membre vous permet d'avoir des options supplémentaires.