Bienvenue, Invité. Merci de vous connecter ou de vous inscrire. Avez-vous oublié d'activer ?

Voir les contributions

Cette section vous permet de consulter les contributions (messages, sujets et fichiers joints) d'un utilisateur. Vous ne pourrez voir que les contributions des zones auxquelles vous avez accès.


Sujets - Kuhn

Pages: [1]
1
Général / Chaine de Toda
« le: Mars 14, 2019, 01:51:48 pm »
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;
}
 

Pages: [1]
anything