Simulation de pendules couples
Description
Voici un petit pack de sources en rapport avec la simulation de pendules couples, je vous le laisse decouvrir par vous meme. Je vous conseille toutefois de pas trop faire attention a la fonciton 'jacobi', il est suffisant de savoir ce qu'elle fait et de savoir l'utiliser (elle diagonalise une matrice et nous permet d'avoir acces aux valeur propres et aux vecteurs propres si je me souvient bien).
J'espere rendre service a tout(e?)s les etudiant(e?)s interesser (ou pas) par ce genre de simulation.
NB: Ce code a etait fait en colaboration avec 'le Grec'.
J'espere rendre service a tout(e?)s les etudiant(e?)s interesser (ou pas) par ce genre de simulation.
Source / Exemple :
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l]; a[i][j]=g-s*(h+g*tau);a[k][l]=h+s*(g-h*tau);
#define C 2 //si jamais on veux remplacer C dans le reste du prog il suffit decrire un autre chiffre que 2
#define W1 1 //idem
double *vecteur(int n);
void ecrire_vecteur(double *v,int n);
double **matrice(int m,int n);
void ecrire_matrice(double **A, int m, int n);
void jacobi(double **a,int n,double *d,double **v1, int *nrot);
void eigsrt(double *d, double **v, int n);
FILE *n_pendules_couples, *fopen(); //fopen est une fonction on doit donc la marquer
main()
{
double *v; // le vecteur
int i,j,k,n,m,nrot;
double **A,*d,**v1,n2;
n_pendules_couples=fopen("n_pendules_couples.dat","w");
//valeurs2=fopen("valeurs2.dat","w");
n=20;
A=matrice(n,n);
d=vecteur(n); //d est un vecteur avec les valeurs propres
v1=matrice(n,n);
for(k=1;k<=n;k++)
{
for(i=1;i<=n;i++)
//init de la matrice... a l'arrache mais ca marche ;-) ...
for(j=1;j<=n;j++)
{
if (i==j)
A[i][j]=-(W1*W1+C);
if (i==j && i==1 && j==1)
{
A[i][j]=-(W1*W1+C);
A[i][j+1]=C;
}
if (i==j && i==n && j==n)
{
A[i][j]=-(W1*W1+C);
A[i][j-1]=C;
}
if (i==j-1)
A[i][j]=C;
if (i==j+1)
A[i][j]=C;
}
//A[1][1]=-(W1*W1+C); // parce que jacobi va me changer la matrice, je reinitialise
//A[1][2]=A[2][1]=C;
//A[2][2]=-(n2*n2+C); // d apres la formule trouvee a part (calcul a la maison)
//ecrire_matrice(A,n,n);
jacobi(A,n,d,v1,&nrot);
eigsrt(d,v1,n);
//printf("valeurs propres: \n"); ecrire_vecteur(d,n);
//printf("vecteurs propres: \n"); ecrire_matrice(v1,n,n);
//printf("\n");
for(i=1;i<=n;i++)
{
for(j=1;j<=n;j++)
fprintf(n_pendules_couples,"%d \t %lg \t",i,v1[i][j]);
fprintf(n_pendules_couples,"\n");
}
}
fclose(n_pendules_couples);
}
// procedure qui permet d' reserver la place pour un vecteur
double *vecteur(int n)
{
double *v;
v=(double*)calloc(n,sizeof(double));
return (v--); // pour que l'indice commence a 1
}
void ecrire_vecteur(double *v,int n) //procedure qui permet d'ecrire un vecteur
{
int i;
for(i=1;i<=n;i++)
printf("%lg \n",v[i]);
}
double **matrice(int m,int n) //declaration de la matrice
{
double**A;
A--; //va changer les indices des lignes
register int i; //pour les utiliser plus rapidement
A=(double**)calloc(m,sizeof(double*)); //on pointe sur les lignes
for(i=1;i<=m;i++)
{
A[i]=(double*)calloc(n,sizeof(double));
A[i]--; //va changer les indices de chaque colonne
}
return A;
}
void ecrire_matrice(double **A, int m,int n) //afficher la matrice
{
int i,j;
for(i=1;i<=n;i++)
{
for(j=1;j<=n;j++)
printf(" %5.3lg ", A[i][j]); // je vais ecrire mes double avec 5 caracteres et 3 chiffres apres la virgule
printf("\n");
}
}
void jacobi(double **a,int n,double *d,double **v, int *nrot) //la matrice va etre modifiee a la sortie, si on veut la garder il faut la sauver
/*Nb: Payes ta perche :D */
{
int j,iq,ip,i;
double tresh,theta,tau,t,sm,s,h,g,c,*b,*z;
b=vecteur(n);
z=vecteur(n);
for (ip=1;ip<=n;ip++){
for (iq=1;iq<=n;iq++) v[ip][iq]=0.0;
v[ip][ip]=1.0;
}
for (ip=1;ip<=n;ip++){ // initialise b et d a la diagonale de a
b[ip]=d[ip]=a[ip][ip]; //ce vecteur va accumuler les termes en tant qu'equations
z[ip]=0.0;
}
*nrot=0;
for (i=1;i<=50;i++) {
sm=0.0;
for (ip=1;ip<=n-1;ip++){ // somme des elements diagonaux
for (iq=ip+1;iq<=n;iq++)
sm += fabs(a[ip][iq]);
}
if (sm == 0.0){
z++;
free(z);
b++;
free(b);
return;
}
if (i < 4)
tresh=0.2*sm/(n*n);
else
tresh=0.0;
for (ip=1;ip<=n-1;ip++){
for (iq=ip+1;iq<=n;iq++){
g=100.0*fabs(a[ip][iq]);
if (i > 4 && (double)(fabs(d[ip])+g) == (double)fabs(d[ip]) && (double)(fabs(d[iq])+g) == (double)fabs(d[iq]))
a[ip][iq]=0.0;
else if (fabs(a[ip][iq]) > tresh) {
h=d[iq]-d[ip];
if ((double)(fabs(h)+g) == (double)fabs(h))
t=(a[ip][iq])/h;
else{
theta=0.5*h/(a[ip][iq]);
t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
if (theta < 0.0) t = -t;
}
c=1.0/sqrt(1+t*t);
s=t*c;
tau=s/(1.0+c);
h=t*a[ip][iq];
z[ip] -= h;
z[iq] += h;
d[ip] -= h;
d[iq] += h;
a[ip][iq]=0.0;
for (j=1;j<=ip-1;j++){
ROTATE(a,j,ip,j,iq)
}
for (j=ip+1;j<=iq-1;j++){
ROTATE(a,ip,j,j,iq)
}
for (j=iq+1;j<=n;j++){
ROTATE(a,ip,j,iq,j)
}
for (j=1;j<=n;j++){
ROTATE(v,j,ip,j,iq)
}
++(*nrot);
}
}
}
for (ip=1;ip<=n;ip++){
b[ip] += z[ip];
d[ip]=b[ip];
z[ip]=0.0;
}
}
printf("too many iterations in routine jacobi \n");
}
//EOD
void eigsrt(double *d, double **v, int n)
{
int k,j,i;
double p;
for (i=1;i<n;i++)
{
p=d[k=i];
for (j=i+1;j<=n;j++)
if (d[j] >= p) p=d[k=j];
if (k != i)
{
d[k]=d[i];
d[i]=p;
for (j=1;j<=n;j++)
{
p=v[j][i];
v[j][i]=v[j][k];
v[j][k]=p;
}
}
}
}Conclusion :
Bon ne m'en veuillez pas trop si la presentation et les commentaires ne sont pas parfaits.NB: Ce code a etait fait en colaboration avec 'le Grec'.
Codes Sources
A voir également
- Simulation de pendules couples
- Simulation keno - Forum Visual Basic
- Simulation schema electrique - Forum VB.NET
- Simulateur de dés ✓ - Forum C++ & C++ .NET
- Simulateur taux de marge ✓ - Forum Javascript / DHTML
- Simulateur de click - Forum Visual Basic
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.