Bonjour , c'est la première fois que je poste sur ce forum .
Voilà , ça fait quelque jours je bosse sur mon programme et il ne fonctionne toujours pas , je commence à en avoir vraiment marre.

Je dois coder en C++ le schema de Runge-Kutta d'ordre 4 pour résoudre le problème dy/dt = f(y) pour 2 ou 3 corps. Mon vecteur y représente les positions et les vitesses des corps dans un plan, y=(x1,y1,x2,y2,x3,y3,vx1,vy1,v x1,vy2,vx3,vy3).

Le problème c'est que lorsque j'étudie la convergence de mes résutats, mon schéma converge d'ordre 1 , ce qui n'est pas normal, mais je ne trouve pas le problème.

Pour étudier la convergence , on simule un astéroide s'approchant de la terre , avec certaine condition initiale , on calcule analytiquement la distance minimale entre l'astéroide et la terre dans un certain intervalle de temps et des pas de temps différents.
Ainsi , on calcule la différence obtenue de chaque distance minimale entre la solution numérique et analytique et on la trace en fonction du pas de temps.

Voici mon schéma de Runge Kutta :



Code:

void PhysEngineRungeKutta::step(){
	
	// y1 est l'étape d'integration avec dt et y2 est l etape d'integration avec dt/2


	
	valarray<double> k1 ;
	valarray<double> k2 ;
	valarray<double> k3 ;
	valarray<double> k4 ;
	bool bon (false);
	valarray<double> y1 ;
	valarray<double> y2;
	valarray<double> r1(2);// terre
	valarray<double> r2 (2);//lune
	valarray<double> r3 (2);//asteroide
	

	
	
	//etape avec dt normale
	k1=dt*acceleration(y);
	k2=dt*acceleration(y+0.5*k1);
	k3=dt*acceleration(y+0.5*k2);
	k4=dt*acceleration(y+k3);
	y=y+(1.0/6.0)*(k1+(2.0*k2)+(2.0*k3)+k4);
}

et acceleration et la fonction f(y) :

Code:

valarray<double> PhysEngine::acceleration (valarray<double> y1)
{
	
	valarray<double> a ;
	valarray<double> r1(2);// terre
	valarray<double> r2 (2);//lune
	valarray<double> r3 (2);//asteroide
	
	
	r1=y1[std::slice(0,2,1)];
	r2=y1[std::slice(2,2,1)];
	r3=y1[std::slice(4,2,1)];
	
	  valarray<double> flag = {1.,1.,1.,1.,0.,0.,1.,1.,1.,1.,0.,0.};

	double G (6.674*pow(10,-11));
	
	a= { y[6] , y[7] , y[8], y[9] , y[10] , y[11] ,
		-G*M_2/pow(distance(r1,r2),2.0)*direction(r1,r2)[0]-    G*M_3/pow(distance(r1,r3),2.0)*direction(r1,r3)[0],
		-G*M_2/pow(distance(r1,r2),2.0)*direction(r1,r2)[1]-    G*M_3/pow(distance(r1,r3),2.0)*direction(r1,r3)[1],
		
		-G*M_1/pow(distance(r2,r1),2.0)*direction(r2,r1)[0]-    G*M_3/pow(distance(r2,r3),2.0)*direction(r2,r3)[0],
		-G*M_1/pow(distance(r2,r1),2.0)*direction(r2,r1)[1]-    G*M_3/pow(distance(r2,r3),2.0)*direction(r2,r3)[1],
		
		-G*M_1/pow(distance(r3,r1),2.0)*direction(r3,r1)[0]-    G*M_2/pow(distance(r3,r2),2.0)*direction(r3,r2)[0],
		-G*M_1/pow(distance(r3,r1),2.0)*direction(r3,r1)[1]-    G*M_2/pow(distance(r3,r2),2.0)*direction(r3,r2)[1]};
	

	
	return a ;
	
}