Chaine de Toda

[Kuhn]

Bonjour, je travaille sur un projet et je suis bloqué c'est pourquoi je viens solliciter votre aide. Je veux représenter le mouvement de trois points qui se déplacent sur un cercle et qui sont en interactions entre eux. Pour ce faire j'ai intégrer numériquement le mouvement des points avec deux méthodes différentes : Runge-Kuta puis Verlet, je code avec C++.
Je demande au code de me fournir un fichier texte puis je représente l'ensemble avec Sfml, j'ouvre une fenêtre dans laquelle crée trois cercles et je les fait tourner en lisant le fichier texte qui va donner l'angle pour la rotation.

Malheureusement ça ne me donne pas du tout le résultat attendu, je ne sais pas si ce sont les méthodes d'intégrations qui sont mauvaises ou bien si j'utilise mal les outils de sfml (ce qui est très possible car j'ai appris à les utiliser depuis deux semaines seulement et justement pour représenter ces points) donc voilà deux semaines que je me casse la tête dessus.
Si vous avez des indications à me proposer le code est ci-dessous :



Intégration

#include <iostream>
#include <fstream>
#include <cmath>
#include <vector>
#include <cstdlib>
#include <math.h>
#define PI 3.14159265
using namespace std;



//Variables globales
double h=0.001;//Pas d'intégration
int pointnumber=3; //Nombre de points dans la chaîne de Toda
double k_phi=0.001; //Force du potentiel
//Variables globales utilisées pour le calcul
vector<double> omega_half(pointnumber);
vector<double> omegapoint_half(pointnumber);




//Fonctions
double F(double omega,double omega_moins,double omega_plus);
void Verlet(vector<double>& omega, vector<double>& omegapoint, vector<double>& omegapointpoint);
void RK2(vector<double>& omega, vector<double>& omegapoint, vector<double>& omegapointpoint);
void update_acceleration(vector<double> omega, vector<double>& omegapointpoint);


//Main
int main()
{
    //Déclaration vecteurs données
    vector<double> omega;
    omega.reserve(pointnumber);
    omega.resize(0);
    vector<double> omegapoint;
    omegapoint.reserve(pointnumber); //donne la capacité du vecteur
    omegapoint.resize(0); //Permet de vider le vecteur
    vector<double> omegapointpoint;
    omegapointpoint.reserve(pointnumber);
    omegapointpoint.resize(0);
   
    omega_half.reserve(pointnumber);
    omegapoint_half.reserve(pointnumber); //donne la capacité du vecteur
   
   
    //Ouverture fichier données
    ofstream valeurs("valeursverlet.txt");
   
   
    //Déclaration conditions initiales
    for (int i=0;i<pointnumber;i++) //Equilibre
    {
        omega.push_back(2*PI/3); //Def d'omega_0
        omegapoint.push_back(0); //Def d'omegapoint_0
    }
    omegapoint[0]=0.01;
    //Remplissage premières lignes fichier
    valeurs<<0<<" "; //Numéro de ligne
    for(int i=0;i<pointnumber;i++)
    {
        omegapointpoint.push_back(0); //Acceleration initiale nulle
        valeurs<<omega[i]<<" "<<omegapoint[i]<<" "<<omegapointpoint[i]<<" ";
    }
    valeurs<<" "<<endl;
    //Définition fenêtre intégration
    double a,b;
    a=0;
    b=1000; //bornes de l'intégrale
    float N=(b-a)/h; //Nombre de points de données
    for (int i=0;i<N;i++)
    {
        RK2(omega, omegapoint, omegapointpoint);
        valeurs<<i+1<<" ";
        for(int j=0;j<pointnumber;j++) //Remplissage fichier
        {
            valeurs<<omega[j]<<" "<<omegapoint[j]<<" "<<omegapointpoint[j]<<" ";
        }
        valeurs<<endl;
    }
    //system("gnuplot -persist plot_erreur.gnu"); //Plot des données
    system("gnuplot -persist portrait_de_phase.gnu");
    return 0;
}

//-----------------------------------------------------------------------

//Fonction reliant position et accélaration dans l'équa diff
double F(double omega,double omega_moins,double omega_plus)
{
    return exp(-(omega-omega_moins))-exp(-(omega_plus-omega));
}


//Permet de calculer les accélérations de tous les points
 void update_acceleration(vector<double> omega, vector<double>& omegapointpoint)
{
    omegapointpoint[0]=(k_phi*F(omega[0],omega[pointnumber-1],omega[1]));
    for (long int i=1;i<pointnumber;i++)
    {
        omegapointpoint.at(i)=(k_phi*F(omega[i],omega[i-1],omega[i+1]));
    }
}


//schéma d'intégration Verlet, update position et accélération, symplectique
void Verlet(vector<double>& omega, vector<double>& omegapoint, vector<double>& omegapointpoint)
{
    for (int i=0;i<pointnumber;i++)
    {
        omegapoint.at(i)+=(h/2)*omegapointpoint.at(i);
    }
    //Schéma Verlet
    for (int i=0; i<pointnumber;i++)
    {
        omega.at(i)=fmod(omega[i]+h*omegapoint[i],2*PI);
    }
    //Permet de calculer les accélérations de tous les points
    update_acceleration(omega,omegapointpoint);
for (int i=0;i<pointnumber;i++)
    {
        omegapoint.at(i)+=(h/2)*omegapointpoint[i];
    }
    return ;
}
void RK2(vector<double>& omega, vector<double>& omegapoint, vector<double>& omegapointpoint)
{
    //Calcul du midpoint
    for (int i=0;i<pointnumber;i++)
    {
        omega_half[i]=fmod(omega[i]+0.5*h*omegapoint[i],2*PI);
        omegapoint_half[i]=omegapoint[i]+0.5*h*omegapointpoint.at(i);
    }
    update_acceleration(omega_half,omegapointpoint);
    for (int i=0;i<pointnumber;i++)
    {
        omegapoint[i]+=h*omegapointpoint.at(i);
    }
    for (int i=0;i<pointnumber;i++)
    {
        omega[i]=fmod(omega[i]+h*omegapoint.at(i),2*PI);
    }
    return;
}

Sfml

#include <SFML/Graphics.hpp>
#include <SFML/Window.hpp>
#include <SFML/System.hpp>
#include <iostream>
#include <string>
#include <fstream>
#include <cmath>
#include <vector>
#include <cstdlib>
#include <math.h>
#include <limits>
#define PI 3.14159265
int main()
{
    //Partie récupération données du fichier
    std::ifstream valeurs("../../projet/valeursverlet.txt");
    if (!valeurs){
        std::cerr<<"Impossible d'ouvrir le fichier"<<std::endl;
        exit(1);
    }
    std::string ligne;
    int pointnumber=3; //Nombre de points dans la chaîne
    int size=0;
    while(valeurs.ignore(std::numeric_limits<int>::max(),'\n')){size++;};
    valeurs.clear();
    valeurs.seekg(0,std::ios::beg); //Revenir au début du fichier
    //Angles
    std::vector<std::vector<double>> omega(size, std::vector<double>(pointnumber));
   
    //Vitesses
    std::vector<std::vector<double>> omegapoint(size, std::vector<double>(pointnumber));
   
    //Accélérations
    std::vector<std::vector<double>> omegapointpoint(size, std::vector<double>(pointnumber));
    std::vector<double> omegatemp(pointnumber);
    std::vector<double> omegapointtemp(pointnumber);
    std::vector<double> omegapointpointtemp(pointnumber);
    int ignore;
        for (int j=0;j<size;j++)
        {
            valeurs>>ignore; //On ignore le numéro de ligne
            for (int i=0;i<pointnumber;i++)
                {
                    valeurs>>omegatemp[i]>>omegapointtemp[i]>>omegapointpointtemp[i];
                    omega[j][i]=omegatemp[i];
                    omegapoint[j][i]=omegapointtemp[i];
                    omegapointpoint[j][i]=omegapointpointtemp[i];
                }
            valeurs.ignore(std::numeric_limits<int>::max(),'\n'); //On passe à la prochaine ligne
        }
   
   
   
    //----------------------------------------------------------------
   
    //Partie vidéo
    sf::ContextSettings settings;
    settings.antialiasingLevel = 8; //permet l'anticrenelage dans la fenetre
   
    float circle_radius=50;
    float fen_size=700;
    sf::RenderWindow window(sf::VideoMode(fen_size, fen_size), "Toda", sf::Style::Default, settings); //cree une fenetre 700/700 pixel? avec un titre ,, settings renvoit au settings precedent pour permetttre l'anticrenelage
    sf::CircleShape circle1(circle_radius); //cree une forme cercle avec pour nom cercle de rayon 50 pixel?
    sf::CircleShape circle2(circle_radius); //cree une forme cercle avec pour nom cercle de rayon 50 pixel?
    sf::CircleShape circle3(circle_radius);
    circle1.setFillColor(sf::Color::Green);
    circle2.setFillColor(sf::Color::Red);
    circle3.setFillColor(sf::Color::Color(255,165,0));//orange
    circle1.setPointCount(1000); // cercle avec 1000 points sur le contour, si on met 3 points on a un triangle
    circle2.setPointCount(1000); // cercle avec 1000 points sur le contour, si on met 3 points on a un triangle
    circle3.setPointCount(1000);
    std::vector<double> x0(pointnumber);
    std::vector<double> y0(pointnumber);
    std::vector<double> theta0(pointnumber);
    for (int i=0;i<pointnumber;i++)
    {
        theta0[i]=fmod(omega[0][i]+(i-1)*2*PI/3,2*PI)*180/PI;
        y0[i]=circle_radius*sin(theta0[i]);
        x0[i]=circle_radius*cos(theta0[i]);
    }
    circle1.setOrigin(-fen_size/2,-fen_size/2);// centre le cercle a l'origine/ setOrigin(0,0) place le cercle dans le coin gauche
    circle2.setOrigin(-fen_size/2,-fen_size/2);
    circle3.setOrigin(-fen_size/2,-fen_size/2);
    sf::Transform t1; //definit une transformation qu'on appelle t1
    sf::Transform t2;
    sf::Transform t3;
    double angle1 = theta0[0] ;
    double angle2 = theta0[1] ;
    double angle3 = theta0[2] ;
    sf::Clock clock;
    sf::Time t0 = sf::seconds(1/60.f);
    int i=0;
    t1.rotate(angle1,350,350);
    t2.rotate(angle2,350,350);
    t3.rotate(angle3,350,350);
    while (window.isOpen() && i<size-1)
    {
        clock.restart();
        while (clock.getElapsedTime()<t0)
        {
            int j=0;
            sf::Event event;
            while (window.pollEvent(event))
            {
                if (event.type == sf::Event::Closed)
                window.close();
            }
            window.clear();
            window.draw(circle1,t1); // dessine le cercle 1 et lui applique t1
            window.draw(circle2,t2); //dessine le cercle et lui applique t2
            window.draw(circle3,t3);
            window.display();
            t1.rotate(angle1*t0.asSeconds(),350,350);
            t2.rotate(angle2*t0.asSeconds(),350,350);
            t3.rotate(angle3*t0.asSeconds(),350,350);
        }
        i++;
        angle1=omegapoint[i][0]*180/PI;
        angle2=omegapoint[i][1]*180/PI;
        angle3=omegapoint[i][2]*180/PI;
       
    }
    return 0;
}
 

[Laurent]

Dire "ça marche pas" et poster tout ton code, dont la majeure partie contient des formules mathématiques, ne nous aidera pas vraiment à t'aider. Sois plus précis quant à ton souci (quel est le résultat actuel ? le résultat attendu ?), réduis ton code à ce qui est vraiment pertinent, etc.

Si tu n'es pas sûr de bien utiliser les classes SFML, fais des tests simples dans un projet à part, sans impliquer les formules mathématiques dont tu n'es pas sûr.

[Kuhn]

Oui vous avez tout à fait raison, c'était juste que je ne savais pas par quoi commencer. Si on regarde le code sfml, posté en second dans mon message

{
    //Partie récupération données du fichier
    std::ifstream valeurs("valeursverlet.txt");
    if (!valeurs){
        std::cout<<"Impossible d'ouvrir le fichier"<<std::endl;
        exit(1);
    }

Je récupère d'abord les valeurs du fichier texte. Le fichier texte se présente comme un tableau, avec en première colonne le numéro de la ligne, la deuxième donne la position (l'angle) du point 1, la troisième donne la vitesse du point 1, la quatrième donne l'accélération du point 1 et etc ... ça donne qqchose comme :
n p1 v1 a1 p2 v2 a2 p3 v3 a3

Ensuite je récupère la position (ainsi que la vitesse et l'accélération) des points que je stock dans des vecteurs séparés :

    for (int j=0;j<size;j++)
        {
            valeurs>>ignore; //On ignore le numéro de ligne
            for (int i=0;i<pointnumber;i++)
                {
                    valeurs>>omegatemp[i]>>omegapointtemp[i]>>omegapointpointtemp[i];
                    omega[j][i]=omegatemp[i];
                    omegapoint[j][i]=omegapointtemp[i];
                    omegapointpoint[j][i]=omegapointpointtemp[i];
                }
            valeurs.ignore(std::numeric_limits<int>::max(),'\n'); //On passe à la prochaine ligne
    for (int i=0;i<pointnumber;i++)
    {
        theta0[i]=fmod(omega[0][i]+(i-1)*2*PI/3,2*PI)*180/PI;
    }

    double angle1 = theta0[0] ;
    double angle2 = theta0[1] ;
    double angle3 = theta0[2] ;
        }

Ensuite je crée une fenêtre vidéo de 700 pixels ainsi que trois cercles qui représentent les trois points (je pense que y a pas de problème là dessus) et je crée des transformations t1 t2 et t3 qui sont des rotate de la classe transform.

 
    circle1.setOrigin(-fen_size/2,-fen_size/2);
    circle2.setOrigin(-fen_size/2,-fen_size/2);
    circle3.setOrigin(-fen_size/2,-fen_size/2);
    sf::Transform t1; //definit une transformation qu'on appelle t1
    sf::Transform t2;
    sf::Transform t3;  
    t1.rotate(angle1,350,350); //fait une rotation de centre de coordonnées (350,350) d'un angle 1
    t2.rotate(angle2,350,350);
    t3.rotate(angle3,350,350);
 

Puis pour gérer le temps entre chaque rotation je crée un variable t0 qui vaut 1/60è de secondes, je mets les rotations et les angles dans une boucle qui va dessiner les cercles et actualiser leurs positions tout en actualisant la valeur des angles.

    while (window.isOpen() && i<size-1)
    {
        clock.restart();
        while (clock.getElapsedTime()<t0)
        {
            int j=0;
            sf::Event event;
            while (window.pollEvent(event))
            {
                if (event.type == sf::Event::Closed)
                window.close();
            }

            t1.rotate(angle1*t0.asSeconds(),350,350);
            t2.rotate(angle2*t0.asSeconds(),350,350);
            t3.rotate(angle3*t0.asSeconds(),350,350);

            window.clear();
            window.draw(circle1,t1); // dessine le cercle 1 et lui applique t1
            window.draw(circle2,t2); //dessine le cercle et lui applique t2
            window.draw(circle3,t3);
            window.display();

        }
        i++;
        angle1=omegapoint[i][0]*180/PI;
        angle2=omegapoint[i][1]*180/PI;
        angle3=omegapoint[i][2]*180/PI;

    }

Le problème c'est que les cercles bougent n'importe comment, ils se passent au travers par exemple alors que c'est pas possible d'après l'intégration qui a l'air assez juste. Tout se passe comme ci le code ne lisait pas le fichier ou plutôt pas comme je le souhaite car le fichier est bien lu (j'ai testé en ajoutant une boucle et ça me donnait bien les 10 premières valeurs du fichiers)

N'hésitez pas à me dire si je dois être plus précis.