////////////////////////////////////////////////////////////////////////////////
// résolution numérique d'une équation de "Schrodinger" 1D 
// -idf/dt=Hf  avec   H=-i(a(x)d/dx+d/dx a(x)=-2ia(x)d/dx-ia'(x)
////////////////////////////////////////////////////////////////////////////////
driver("X11");  
////////////////////////////////////////////////////////////////////////////////
// discrétisation espace 
///////////////////////////////////////////////////////////////////////////////
X=30;// intervalle [-X,X]
dx=0.1;// pas en espace 
K=2*(X/dx)+1;//nombre de points
x=[-X:dx:X]';// vecteur contenant les abcisses
dt=0.002;//pas de temps
///////////////////////////////////////////////////////////////////////////////
// le potentiel v(x) et la donnée initiale u0    
//////////////////////////////////////////////////////////////////////////////
u=%e.^(-(x+5).^2);//une gaussienne
//deff('y=a(x)','y=(1+x.^2)./(2+x.^2)');//le  potentiel
deff('y=a(x)','y=1-0.85*exp(-x.^2)');//le  potentiel
pot=a(x);
A=(a(x));
da=([0;A(1:K-1)]-[A(2:K);0])/(2*dx);//dérivé de a(x)
//////////////////////////////////////////////////////////////////////////////
//lancement affichage
//////////////////////////////////////////////////////////////////////////////
Y=1.5;rect=[-X,-0.5,X,Y];//fenêtre graphique
titlepage(["résolution numérique d EDP";" ";"équation de Schrodinger ";" ";"-idu/dt=Hu, H=a(x).D+D.a(x)"])
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=-Hu@a(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:4000
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*A)/dx).*(rightu-leftu);//+(dt)*(da.*u);
//u(1)=0;u(lon)=0;//condition Dirichlet
u(1)=u(2);u(lon)=u(lon-1);//condition neuman
v=u./A;
norme=sqrt(abs(v));
wave=%e.^(-(x+5-4*i*dt*a(x)).^2);
// affichage du résultat 
xclea(-X,Y,2*X,2*Y)
plot2d([x x],[norme pot], [2 3],"111","-idu/dt=Hu@a(x)@u0(x+ta)",rect)
xstring(10,-0.25,["t=" string(i*dt)])
xset("wshow")
end
xset("pixmap",0)