////////////////////////////////////////////////////////////////////////////////
// résolution numérique d'une équation de Schrodinger 1D 
// -idf/dt=-d²f/dx²+v(x)*f
// phi=phi(t,x)    phi=exp(i(x*xi-t*xi^2))*u    phi|[t=0]=exp(ix*xi-t*xi^2)*u0
////////////////////////////////////////////////////////////////////////////////
driver("X11");  
////////////////////////////////////////////////////////////////////////////////
// discrétisation espace 
///////////////////////////////////////////////////////////////////////////////
X=20;// intervalle [-X,X]
dx=0.08;// pas en espace 
K=2*(X/dx)+1;//nombre de points
x=[-X:dx:X]';// vecteur contenant les abcisses
dt=0.00005;//pas de temps
///////////////////////////////////////////////////////////////////////////////
// le potentiel v(x) et la donnée initiale u0    
//////////////////////////////////////////////////////////////////////////////
u=%e.^(-(x+5).^2);//une gaussienne
deff('y=f(x)','y=bool2s((x<2)&(x>0))');//le support du potentiel
pot=f(x); V=1000*pot;// l'amplitude du potentiel
xi=50;//vitesse
//////////////////////////////////////////////////////////////////////////////
//lancement affichage
//////////////////////////////////////////////////////////////////////////////
Y=2;rect=[-X,-0.5,X,Y];//fenêtre graphique
titlepage(["résolution numérique d EDP";" ";"équation de Schrodinger 1D";" ";"-idu/dt=-d²u/dx²+v(x)*u"])
xtitle("(clicquer dans la fenêtre pour lancer l animation)")
[c_i,c_x,c_y,c_w]=xclick();//attendre un clicque pour passer aux graphes
xbasc();
xset("pixmap",1) 
xtitle("équation de Schrodinger 1D (clicquer pour lancer)","x","amplitude")
norme=abs(u);
lon=length(u);
plot2d([x x],[norme pot], [2 3],"111","-idu/dt=-d²u/dx²+v(x)u@v(x)",rect)
xset("wshow")
[c_i,c_x,c_y,c_w]=xclick();//attendre un click pour démarer
/////////////////////////////////////////////////////////////////////////////////
//calcul récursif de l'évolution 
/////////////////////////////////////////////////////////////////////////////////
for i=1:2300
leftu=[u(2:K); 0];//u(k+1)
rightu=[0;u(1:K-1)];//u(k-1)
// calcul de u au pas de temps suivant 
u=u+(2*(dt*xi)/dx)*(rightu-leftu)+((%i*dt)/(dx^2))*(2*u-(leftu+rightu))-(%i*dt)*(V.*u);
//u(1)=0;u(lon)=0;//condition Dirichlet
u(1)=u(2);u(lon)=u(lon-1);//condition neuman
norme=abs(u);
// affichage du résultat 
xclea(-X,Y,2*X,2*Y)
plot2d([x x],[norme pot], [2 3],"111","-idu/dt=-d²u/dx²+v(x)u@v(x)",rect)
xstring(10,-0.25,["t=" string(i*dt)])
xset("wshow")
end
xset("pixmap",0)