1. Accueil
  2. C / C++ / C++.NET
  3. Sources C / C++ / C++.NET
  4. Sources Maths & Algorithmes
  5. Solution graphique approchée au problème des n corps

Solution graphique approchée au problème des n corps

Soyez le premier à donner votre avis sur cette source.

Vue 8 740 fois - Téléchargée 277 fois

Télécharger le projet
Commenter

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.

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;

			break;
		case 2:


    			
  • (Uint16 *)p = (Uint16) color;

			break;
		case 3:
			if(SDL_BYTEORDER == SDL_BIG_ENDIAN)
			{


    				
  • (Uint16 *)p = ((color >> 8) & 0xff00) | ((color >> 8) & 0xff);

    • 				
  • (p+2) = color & 0xff;

			}
			else
			{


    				
  • (Uint16 *)p = color & 0xffff;

    • 				
  • (p+2) = (color >> 16) & 0xff;

			}
			break;
		case 4:


    			
  • (Uint32 *)p = color;

			break;
	}
}

/* Anciennes positions des corps dans le systeme de coordonnees de l'ecran */
Uint16 old_pos_x[N];
Uint16 old_pos_y[N];

void animer(corps *tab, int n, double dt, SDL_Surface *screen)
{
	int i;
	int need_update = 0;
	while(1)
	{
		maj(tab,n,dt);

		for(i=0;i<n;i++)
		{	
			Uint16 new_x = (Uint16)tab[i].pos.x;
			Uint16 new_y = (Uint16)tab[i].pos.y;
			if((new_x <= 0) || (new_x >= XRES)
			 ||(new_y <= 0) || (new_y >= YRES))
			{
			 	printf("Un corps a quitté l'écran. Fin du programme\n");
				return;
			}
			/* On ne change l'affichage que si le corps a effectivement bouge a l'ecran */
			if((new_x != old_pos_x[i]) || (new_y != old_pos_y[i]))
			{
				old_pos_x[i] = new_x;
				old_pos_y[i] = new_y;
				SDL_LockSurface(screen);
				PutPixel(screen,new_x, new_y, color[i % NB_COLOR]);
				SDL_UnlockSurface(screen);
				need_update++;
			}
		}
		/* Pour pallier a la lenteur de l'affichage on ne le rafraichit que toutes les 10 modifications */
		if(need_update >= 10)
		{
			need_update = 0;
			SDL_Flip(screen);
		}	
	}
}

int main(int argc, char **argv)
{
	SDL_Surface *screen;
	corps tab[]= {	{1000000000 ,{512,384},	vecteur_nul	,vecteur_nul},		/* le "soleil" */
			{1000000    ,{350,384},	{0,0.01}	,vecteur_nul},		/* les planetes */
			{3000000    ,{512,450},	{0.02,0} 	,vecteur_nul},
			{1000500    ,{300,300},	{0.007,-0.001}	,vecteur_nul},
			{2000000    ,{650,450},	{-0.003,0.003}	,vecteur_nul}};

	if(SDL_Init(SDL_INIT_VIDEO) == -1)
	{
		printf("Initialisation de SDL impossible : %s\n",SDL_GetError());
		return -1;
	}
	screen = SDL_SetVideoMode(XRES,YRES,24,SDL_VIDEO_FLAGS);
	SDL_FillRect(screen, NULL, SDL_MapRGB(screen->format, 0x00,0x00,0x00));	// écran noir
	SDL_ShowCursor(SDL_DISABLE);
	SDL_Flip(screen);
	
	animer(tab,N,1,screen);
	
	SDL_Quit();	
	return 0;
}

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.
Télécharger le projet

Codes Sources

A voir également

Ajouter un commentaire Commentaires
cs_JCDjcd Messages postés 1138 Date d'inscription mardi 10 juin 2003 Statut Membre Dernière intervention 25 janvier 2009 4
22 mai 2005 à 13:45
Programme amusant, moi j'en ai fait un pareil, mais je ne l'ai pas encore mis, je dois encore l'ameliorer, il faut que j'implemente du Runge-Kutta, car ici du fait du Euler, or le probleme est numeriquement instable.

Note historique : Poincaré a demontrer qu'il n'existait pas de solution générale, donc pas besion de chercher ...

Je trouve le code bien fait (partitionnement des fonctions) et bien presenter.
MetalDwarf Messages postés 241 Date d'inscription mardi 29 octobre 2002 Statut Membre Dernière intervention 23 janvier 2006
22 mai 2005 à 19:16
En fait je l'ai fait en Caml au début, tout du moins le début, c'était une partie d'un sujet d'informatique de Polytechnique il y a quelques annees (il etait question d'un approximation a base de quadtree permettant de passer le probleme en O(n.log(n))).
C'est vrai que le probleme a des chances d etre chaotique, c est pour ca que le "dt" est petit (0.1s) pour eviter les erreurs.
Faudra que je planche sur Runge-Kitta un jour, Euler c'est un peu limite en effet.
Jarod1980 Messages postés 273 Date d'inscription samedi 5 juillet 2003 Statut Membre Dernière intervention 31 mars 2015 2
23 mai 2005 à 17:34
Très bon code et très bien commenté! en ce qui concerne Runge Kutta il y en a sur le site dont un de moi je crois que JCDjcd en a publié un également.
Des systèmes d'une telle complexité sont chaotiques. Pour N>3 ils n'existent pas de solution générale.
cs_Kirua Messages postés 3006 Date d'inscription dimanche 14 avril 2002 Statut Membre Dernière intervention 31 décembre 2008
23 mai 2005 à 17:42
Euhm, tu l'as sans doute fait et tu n'as pas précisé la constante du n² comme c'est l''habitude, mais tu peux passer ton algo en O(n²/2) en considérant les couples de corps, comme les forces en jeux entre ces deux corps sont simplement opposées. Ceci dit, je dis pê qq ch que tu as fait, désolé ;).
MetalDwarf Messages postés 241 Date d'inscription mardi 29 octobre 2002 Statut Membre Dernière intervention 23 janvier 2006
23 mai 2005 à 18:40
Jarod1980 > Oui j'ai deja vu des choses la dessus sur le site, mais je regarderais ca quand j aurais le temps (apres les concours, et oui math spe c'est pas drole tous les jours)

Kirua > Tres bonne idee, d'ailleurs c'est fait! (depuis 10 minutes ;). Cependant la plus grande part du temps processeur est "mangee" par SDL, si vous savez comment changer ca... Je ne connais absolument rien a SDL.
S'inscrire maintenant

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.

Du même auteur (MetalDwarf)