Solution graphique approchée au problème des n corps
Description
Le titre est assez explicite je pense, sauf pour ceux qui ne savent pas ce que c'est... Le problème des N corps est un problème de physique qui consiste à calculer les trajectoires de corps soumis a l'attraction gravitationnelle. Ce problème n'a pas de solution "analytique", c'est à dire exacte pour N>2 (du moins personne ne l'a encore trouvée), d'ou l'importance de le résoudre de façon approchée. La solution présentée ici est en 2D, et fait intervenir les relations de base de la mécanique (le PFD, principe fondamental de la dynamique).
Le code n'est pas optimisé au maximum, en particulier certains calculs pourraient etre evites (comme le calcul de la force alors que celui de l ecceleration peut etre fait directement), mais est "collé" aux formules physique correspondantes.
L'affichage est géré par SDL, c'est mon premier code en SDL, et d'ailleurs la fonction d'affichage d'un pixel est reprise d'un Linux Magazine, mais l'interet n'est pas la. Ce programme a été testé sous Linux avec SDL 1.2, et peut etre assez lent sur un PC peu puissant (dans ce cas augmenter la valeur de dt, qui est la precision de la simulation). Enfin, L'algorithme utilisé est FONDAMENTALEMENT un algorithme en O(N²) par rapport au nombre de corps, donc ca peut devenir lent pour de tres nombreux corps.
Je crois que le code n'est pas trop dur à comprendre, n'hesitez pas a formuler des commentaires.
Pour rajouter des corps, il suffit de changer la valeur de N, et de les definir dans main(). Ce n est pas tres evolue, mais c est un petit code fait en 2h.
Le systeme de corps present dans la source montre que les trajectoires des "planetes" autour d'un "soleil" sont des ellipses, mais seulement en premiere approximation (si on considere la masse du soleil infinie par rapport a celle des planetes et qu'on neglige l'attraction entre les planetes devant celle exercee par le soleil). C'est pour cela que les trajectoires ne sont pas périodiques.
Le code n'est pas optimisé au maximum, en particulier certains calculs pourraient etre evites (comme le calcul de la force alors que celui de l ecceleration peut etre fait directement), mais est "collé" aux formules physique correspondantes.
L'affichage est géré par SDL, c'est mon premier code en SDL, et d'ailleurs la fonction d'affichage d'un pixel est reprise d'un Linux Magazine, mais l'interet n'est pas la. Ce programme a été testé sous Linux avec SDL 1.2, et peut etre assez lent sur un PC peu puissant (dans ce cas augmenter la valeur de dt, qui est la precision de la simulation). Enfin, L'algorithme utilisé est FONDAMENTALEMENT un algorithme en O(N²) par rapport au nombre de corps, donc ca peut devenir lent pour de tres nombreux corps.
Source / Exemple :
/* Solution graphique approchée du problème des N corps */
/* Metaldwarf 2005 */
/* Testé sous Linux avec la libSDL 1.2 */
/* gcc ncorps.c -o ncorps -lm -Wall `sdl-config --cflags --libs` */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <SDL/SDL.h>
#define SDL_VIDEO_FLAGS (SDL_HWSURFACE | SDL_ANYFORMAT)
/* Le nombre de corps */
#define N 5
/************************ Partie Algorithmique ************************/
/* La constante universelle de gravitation */
#define G 6.67e-11
/* vecteur du plan */
typedef struct _vect
{
double x;
double y;
} vect;
vect vecteur_nul = {0,0};
typedef struct _corps
{
double m; /* masse */
vect pos; /* position */
vect v; /* vitesse */
vect a; /* acceleration */
} corps;
/* norme d'un vecteur */
double norme2(vect u)
{
return sqrt((u.x * u.x) + (u.y * u.y));
}
/* res = u1 + u2 */
vect vect_add(vect u1, vect u2)
{
vect res;
res.x = u1.x + u2.x;
res.y = u1.y + u2.y;
return res;
}
/* res = u1 - u2 */
vect vect_sub(vect u1, vect u2)
{
vect res;
res.x = u1.x - u2.x;
res.y = u1.y - u2.y;
return res;
}
/* res = a.u (a réel) */
vect mul_scal(double a, vect u)
{
vect res;
res.x = a * u.x;
res.y = a * u.y;
return res;
}
/* force gravitationnelle exercée par p1 sur p2 */
vect force_grav(corps p1, corps p2)
{
double d; /* distance entre les deux corps */
vect u; /* vecteur unitaire dirigé de p2 vers p1 */
d = norme2(vect_sub(p1.pos, p2.pos));
u = mul_scal((1 / d),vect_sub(p1.pos,p2.pos));
return mul_scal(G * (p1.m * p2.m) * (1/(d*d)), u);
}
/* application du principe fondamental de la dynamique au point p, d'index i dans un tableau de taille p */
/* modifie son accélération */
void pfd(corps *tab, int i, int n)
{
int k;
vect f = vecteur_nul; /* vecteur force */
corps p = tab[i];
for(k=0;k<n;k++)
{
if(k != i)
f = vect_add(f, force_grav(tab[k], p));
}
tab[i].a = mul_scal((1 / p.m),f);
}
/* optimisation : F(1->2) = -F(2->1) donc on peut éviter la moitie des calculs */
void pfd2(corps *tab, int i, int n)
{
int k;
corps p = tab[i];
vect v = vecteur_nul;
/* premier corps : on reinitialise l'acceleration de tous les corps */
if(i==0)
for(k=0;k<n;k++)
tab[k].a = vecteur_nul;
/* Dernier corps : toutes les contributions ont deja ete calculee */
if(i==(n-1)) return;
for(k=i+1;k<n;k++)
{
v = force_grav(tab[k],p);
tab[i].a = vect_add(tab[i].a,mul_scal((1/p.m),v));
tab[k].a = vect_add(tab[k].a,mul_scal((-1/tab[k].m),v));
}
}
/* on a : dv/dt = a et dpos/dt = v
donc : dv = a * dt et dpos = v * dt */
/* modifie la position d'un corps (on considère que l'accélération a déja été calculée */
void change_pos(corps *tab, int i, double dt)
{
tab[i].v = vect_add(tab[i].v , mul_scal(dt,tab[i].a));
tab[i].pos = vect_add(tab[i].pos , mul_scal(dt,tab[i].v));
}
/* met à jour le tableau de corps en avancant de dt dans le temps */
/* NB : Les deux boucles doivent etre laisees comme ceci sinon l'acceleration des corps n'est pas calculee correctement */
/* En effet certains corps on deja bouges avant de calculer l'acceleration des suivants */
void maj(corps *tab, int n, double dt)
{
int i;
for(i=0;i<n;i++) /* on change les accélérations des corps */
{
pfd2(tab,i,n);
}
for(i=0;i<n;i++) /* puis on met à jour les positions */
{
change_pos(tab,i,dt);
}
}
/************************ Partie Graphique ************************/
#define NB_COLOR 5
#define XRES 1024
#define YRES 768
Uint32 color[NB_COLOR] = {0x00ff0000, 0x0000ff00, 0x000000ff, 0x00ff00ff, 0x00ffff00};
/* repris de Linux Magazine 65 (p.80) */
void PutPixel(SDL_Surface *surface, Uint16 x, Uint16 y, Uint32 color)
{
Uint8 bpp = surface->format->BytesPerPixel;
Uint8 *p = ((Uint8 *) surface->pixels) + y * surface->pitch + x * bpp;
switch(bpp)
{
case 1:
- p = (Uint8) color;
- (Uint16 *)p = (Uint16) color;
- (Uint16 *)p = ((color >> 8) & 0xff00) | ((color >> 8) & 0xff);
- (p+2) = color & 0xff;
- (Uint16 *)p = color & 0xffff;
- (p+2) = (color >> 16) & 0xff;
- (Uint32 *)p = color;
Conclusion :
Je crois que le code n'est pas trop dur à comprendre, n'hesitez pas a formuler des commentaires.
Pour rajouter des corps, il suffit de changer la valeur de N, et de les definir dans main(). Ce n est pas tres evolue, mais c est un petit code fait en 2h.
Le systeme de corps present dans la source montre que les trajectoires des "planetes" autour d'un "soleil" sont des ellipses, mais seulement en premiere approximation (si on considere la masse du soleil infinie par rapport a celle des planetes et qu'on neglige l'attraction entre les planetes devant celle exercee par le soleil). C'est pour cela que les trajectoires ne sont pas périodiques.
Codes Sources
A voir également
- Solution graphique approchée au problème des n corps
- Pedantix solution - Forum C
- Agence interface graphique ✓ - Forum Java
- Logiciel pour créer interface graphique - Forum Java
- Solution pour Pilote Isam Introuvable - Forum Visual Basic
- Tracer graphique en ligne - Forum C / C++ / C++.NET
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.