Simulation d'orbite avec RK4

[dridri85]

Bonjour,
je développe en C un moteur de jeu entièrement de zéro, je souhaite n'utiliser aucune librairie existante.
Du coup j'en arrive à la physique, j'utilise la méthode Runge-Kutta d'ordre 4 (RK4) et elle me donne de très bon résultat pour beaucoup d'utilisations.
Mais j'ai bien envie de faire une démo spatiale, et j'aimerai bien pouvoir intégrer la gravitation newtonienne à mon RK4 plutôt que de faire de la simple et fausse rotation car ça me permettrait par exemple de gérer des collisions entre objets en orbite.

Le problème : la simulation d'orbite ne fonctionne pas très bien, pour le moment je teste avec un satellite seul, il tourne bien autour de ma planète mais son orbite s'agrandit peu à peu et il fini par partir dans l'espace. Bien entendu j'utilise des types 'double' dans mon code pour éviter des erreurs d'arrondi.
Je ne sais pas vraiment si je dois mettre mon code ici puisque c'est tout langages confondus. L'intégrateur RK est tout ce qu'il y a des plus classique, voilà ce qu'il donne en pseudo-code :

Code :

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
deriv Evaluate0(local state)
{
	deriv output
	output.force = state.force
	output.velocity = state.velocity
	return output
}
 
deriv Evaluate(local state, t, dt, d)
{
	state.position = state.position + d.velocity * dt
	state.momentum = state.momentum + d.force * dt
	state.velocity = state.momentum * (1.0 / state->mass)
 
	deriv output
	output.force = state.force
	output.velocity = state.velocity
	return output
}
 
void Integrate(ref state)
{
	a = Evaluate0(*state)
	b = Evaluate(*state, t, dt*0.5f, a)
	c = Evaluate(*state, t, dt*0.5f, b)
	d = Evaluate(*state, t, dt, c)
 
	state->position = state->position + (a.velocity + (b.velocity + c.velocity) * 2.0 + d.velocity) * 1.0f/6.0 * dt
	state->momentum = state->momentum + (a.force + (b.force + c.force) * 2.0 + d.force) * 1.0f/6.0 * dt
 
	state->velocity = state->momentum * (1.0 / state->mass)
}
 
void ApplyForce(ref state, vec)
{
	state->force = state->force + vec
}
 
void ApplyGravity(ref s1, ref s2)
{
	const G = 6.67234E-11
	d = dist(s1->position, s2->position)
	P = G * (s1->mass * s2->mass) / (d * d)
 
	force = P * (s2->position - s1->position) / d
	ApplyForce(s1, force)
	ApplyForce(s2, -force)
}

Voilà comment je calcule la vitesse initiale (d'après Wikipédia : http://fr.wikipedia.org/wiki/Vitesse_orbitale) :

Code :

1
2
satellite.momentum.y = sqrt(G * earth.mass / R) * sat.mass
// sqrt étant la racine carré, G la constante gravitationnelle, et R le rayon orbitale du satellite

Par exemple pour un rayon orbital de 20 000 km, j'ai une vitesse initiale v0 de 4464.18 m/s

Je suppose être dans le vrai si je dis que c'est un problème d'arrondi au niveau des calculs ? Comme faire pour éviter ça sans perdre en performance, y a-t-il des moyens de simplifier les calculs ?
Voilà, si quelqu'un aurait une idée ou un lien magique à propos de ce genre de chose, je n'ai rien trouvé sur google à propos du newtonien avec RK4..

Merci

[pyros]

Salut,

Je ne pense pas que ça vienne de la précision des float/double.

Pour avoir un ordre de grandeur, un double possède 52 chiffre significatif en base 2, ça fait à peut près 16 chiffres significatifs en base 10. Ca fait que sur un système planétaire de dimension environ 1 000 000km (3 fois distance terre/lune), tu as une précision supérieur au micro mètre .

Je pense plus à un problème de formule, de dt trop grand, ou de l'ordre dans lequel tu calcul accélération, vitesse, position. Il y a une ordre à respecter pour que ça marche bien mais je ne me souviens jamais lequel .