Bonjour,

Je me suis depuis peu lancé le petit défi d'arriver à simuler et visualiser l'Équation de Schrodinger 1D stationnaire avec comme librairie d'affichage SFML.

Seulement, j'ai deux problème: quand je met un potentiel escalier, ça dégénère. Et si j'en met pas (ou bien potentiel = 0 partout) bah j'ai ça diverge aussi mais vers 0, et ceci même avec des conditions de Dirichlet (psy(0) = psy(width) = 0).

De plus, problème prépondérant selon moi, ma densité de probabilité n'avance pas vers les x positifs ! Je sais bien qu'il n'y a pas de dépendance du temps puisqu'on prend le module de l'exponentielle complexe qui dépend du temps pour afficher la densité de probabilité. Mais comment le faire marcher autrement ? Je ne sais pas.

Donc voilà ça fait bien 7h que j'essaie de comprendre, si jamais une âme charitable et intéressée par ce mini projet souhaite m'éclairer, c'est avec plaisir !

Voici mon code:

La fonction plot_densité_proba pour plot la densité de probabilité. Si j'utilise pas cette fonction pour plot la courbe, j'utilise "plot_real" qui permet de plot uniquement la partie réelle de la fonction solution de l'équation.

#include <iostream>
#include <SFML/Graphics.hpp>
#include <cmath>
#include <complex>
using namespace std;
using namespace sf;
# define PI 3.14159265358979323846  /* pi */

complex<double> I(0.0, 1.0); //imaginary unit

const int height = 900;
const int width = 900;

const int potential_value = 20;
const double Energie = 22;

const double dx = 1e-2; // 1/400 m
const double dt = 1e-10;
const double mass = 1e-2; //kg
const double hbar = 1e-5; //J.s

//juste pour afficher la densité de probabilité, ne pas prêter attention
auto plot_densite_proba(array<complex<double>, width> vecteur_a_afficher) {
    VertexArray curve(PrimitiveType::LineStrip, std::ceil((width - 0) / 1.f));

    for (auto x = 0; x < width; x += 1) {
        curve.append(Vertex(Vector2f(x, abs(vecteur_a_afficher[x]*conj(vecteur_a_afficher[x])))));
    }

    return curve;
}

//juste pour afficher la partie réelle, ne pas prêter attention
auto plot_real(array<complex<double>, width> vecteur_a_afficher) {
    VertexArray curve(PrimitiveType::LineStrip, std::ceil((width - 0) / 1.f));

    for (auto x = 0; x < width; x += 1) {
        curve.append(Vertex(Vector2f(x, real(vecteur_a_afficher[x])*1000)));
    }

    return curve;
}

//Définition du potentiel
double potentiel(int position) {

    if(position <= width/2 + width/5) {
        return (double)potential_value*exp(-pow((float)(position-width/2-width/5)/50.f,2));
    }
    else {
        return potential_value;
    }

   return 0;

}

int main() {
    //====== PARTIE AFFICHAGE
    RenderWindow window(VideoMode(width + 100, height + 100), "Schrodinger's equation in 1D");
    window.setFramerateLimit(60.f);
    View view; view.setCenter(sf::Vector2f(width/2, height/2)); view.setSize(sf::Vector2f(3*width/2, -2*height)); window.setView(view);
    Transform transform; transform.translate(0, height/2);

    VertexArray densite_proba;
    VertexArray potentiel_draw(PrimitiveType::LineStrip, std::ceil((width - 0) / 1.f));
    for (auto x = 0; x < width; x += 1) {
        Vertex vertexx = Vertex(Vector2f(x, potentiel(x)*potentiel(x)));
        vertexx.color = Color::Blue;
        potentiel_draw.append(vertexx);
    }
    //==== FIN DE LA PARTIE AFFICHAGE

    //Mon psi
    array<complex<double>, width> psy = {0.f};
    
    //Conditions initiales
    for(int x = 0; x < width; x++) {
        psy[x] = Energie*exp(-2*(double)(x - width/5)*(x - width/5)/10000.f);
        psy[x] *= exp(I * (double)x);
    }

    //Boucle
    while (window.isOpen()) { Event event; while (window.pollEvent(event)) { if (event.type == Event::Closed) window.close(); }
        densite_proba = plot_densite_proba(psy);
        //densite_proba = plot_real(psy);

        array<complex<double>, width> n_psy;
        for(int x = 1; x < width - 2; x++) {
            //Équation de schrodinger
            n_psy[x] = psy[x] - I*dt/hbar*(potentiel(x)*psy[x] - hbar*hbar/2/mass/dx/dx * (psy[x + 1] - (double)2*psy[x] + psy[x - 1]));
        }
        n_psy[0] = psy[1]; // Condition périodique à gauche
        n_psy[width - 1] = psy[width - 2]; // Condition périodique à droite

        psy = n_psy;


        //==== Affichage
        window.clear(); window.draw(densite_proba, transform); window.draw(potentiel_draw, transform); window.display(); } return 0;
}

Voilà, des bisous ! Et merci à quiconque répondra

-
Edité par MathiasLoliter 29 septembre 2024 à 9:48:47

Avez-vous fait une analyse de la stabilité numérique de votre algorithme ?

(juste avec des doubles "standard" ? sans bibliothèque de calcul scientifique ?)

Bonjour,

Au niveau de la stabilité, j'ai vérifié que h.dt/2.m.dx² <= 1 ou même que dt/dx² était assez petit. On peut voir dans mon programme que c'est le cas, et même avec une valeur de h.dt/2.m.dx² très petite (1e-5,-6 etc), j'avais tout de même une divergence (et accessoirement l'onde qui n'avancait pas vers les x positifs, qui est aussi un gros problème).

D'ailleurs oui je n'ai pas mentionné ce problème mais quelque chose qui me chagrine est que l'onde n'avance pas. J'ai essayé de mettre une condition initiale sur la dérivée, ou bien de rajouter un terme en exp(ikx) mais ça n'avait pas l'air de plus marcher...

EDIT FINAL - Ok alors en fait avec un peu de recherche, il s'avère qu'on ne peut pas résoudre ce problème avec la méthode d'Euler de façon stable.

J'ai réussi à le simuler, mais il m'a fallut utiliser une autre formule de discrétisation, qui prend en compte une normalisation (pour ne pas avoir de divergence): https://ergodic.ugr.es/cphys/lecciones/SCHROEDINGER/ajp.pdf

-
Edité par MathiasLoliter 29 septembre 2024 à 10:54:04