x_message(['Integrazione problema n.3  con';...
           'schemi alle differenze finite.';' ';...
           ' '; ...
           'Ver 2.1  (11/5/2005)'])



cont0="Y";
while cont0=="Y",


txt=['lambda';'mu    '];
sig=x_mdialog('Parametri equazione campione',txt,['1.0';'0.1']);
lambda=evstr(sig(1));
mu=evstr(sig(2));


N=101;
dt=0.499e-3;
itermax=1000;
Resmin=1e-7;


cambia=2;
while cambia==2,
sN=msprintf('%i', N);
sdt=msprintf('%f', dt);
sitermax=msprintf('%i', itermax);
sResmin=msprintf('%e', Resmin);
txt=['N. nodi';'dt    ';'N. max iter';'Res. min'];
//sig=x_mdialog('Parametri  integrazione ',txt,['101';'0.499e-3';'1000';'1e-7']);
sig=x_mdialog('Parametri  integrazione ',txt,[ sN ; sdt ; sitermax ; sResmin]);
N=evstr(sig(1));
dt=evstr(sig(2));
itermax=evstr(sig(3));
Resmin=evstr(sig(4));

//N=101;
dx=1/(N-1);
c=lambda*dt/dx;
s=mu*dt/dx^2;


tschema=['FTCS (esplicito,  I ord. tempo, II ord. spazio)';...
         'BTCS (implicito,  I ord. tempo, II ord. spazio)';...
         'CN   (implicito, II ord. tempo, II ord. spazio)'];
nschema=x_choose(tschema,['Scegli lo schema alle differenze finite']);


testo(1)=strcat( ['------------------------------------------------------------']);
testo(2) =strcat( ['Riassunto dei parametri equazione ed  integrazione:']);
testo(3) =strcat( ['Parametri equazione ----->  lambda:' string(msprintf('%2.5f',lambda   ))']);
testo(4) =strcat( ['                                mu:' string(msprintf('%2.5f',mu       ))']);
testo(5) =strcat( ['                         mu/lambda:' string(msprintf('%2.5f',mu/lambda))']);
testo(6) =strcat( ['Parametri integrazione --> N. nodi:' string(msprintf('%2.5f',N   ))']);
testo(7) =strcat( ['                                dx:' string(msprintf('%2.5f',dx  ))']);
testo(8) =strcat( ['                                dt:' string(msprintf('%2.5f',dt  ))']);
testo(9) =strcat( ['                                 c:' string(msprintf('%2.5f',c   ))']);
testo(10) =strcat( ['                                 s:' string(msprintf('%2.5f',s   ))']);
testo(11)=strcat( ['                       N. iter max:' string(msprintf('%2.5f',itermax))']);
testo(12)=strcat( ['                         Res. min.:' string(msprintf('%2.5f',Resmin))']);
testo(13)=strcat( ['            Schema di integrazione:' string(msprintf('%s',tschema(nschema) ))']);
testo(14)=strcat( ['------------------------------------------------------------']);
testo(15)=strcat( ['Vuoi cambiare i parametri?']);

cambia= x_message(testo,['No' 'Yes']);
end,




X=[0 : dx : 1.]';
gn=X;
gnp1=zeros(1,N)';

iter=2000;


if nschema==1 then
  //FTCS
  A=zeros(N,N);
  A(1,1)=1.;
  for i=2:N-1
    A(i,i+1)=s-c/2;
    A(i,i)=1-2*s;
    A(i,i-1)=(c/2+s);
  end
  A(N,N)=1.;
  L=A;
elseif nschema==2 then
  //BTCS
  A=zeros(N,N);
  A(1,1)=1.;
  for i=2:N-1
    A(i,i+1)=c/2-s;
    A(i,i)=1+2*s;
    A(i,i-1)=-(c/2+s);
  end
  A(N,N)=1.;
  L=inv(A);
else
  //CN
  A=zeros(N,N);
  B=zeros(N,N);
  A(1,1)=1.;
  B(1,1)=1.;
  for i=2:N-1
    A(i,i+1)=1/2*(c/2-s);
    B(i,i+1)=1/2*(s-c/2);
    A(i,i)=1+s;
    B(i,i)=1-s;
    A(i,i-1)=-1/2*(c/2+s);
    B(i,i-1)=1/2*(c/2+s);
  end
  A(N,N)=1.;
  B(N,N)=1.;
  L=inv(A)*B;
end,




xset("window",0)
xbasc();
xtitle('Evoluzione della soluzione ', 'x', 'g');
xset("window",1)
xbasc();
xtitle('Storia del residuo ', 'iter.', 'Res.');
xset("window",0)
xselect();
kp=xget("pixmap");xset("pixmap",1);
xset("pixmap",kp);
kp=xget("pixmap");xset("pixmap",1);
xselect();plot2d(X,gn,1,rect=[0.,0.,1.,1.]);

xset("mark",-1,2);
t=0.;
res=1e+39;
nn=1;
gnp1max=0.;
t=0.;
while nn<=itermax & res > Resmin, 
  gnp1=L*gn;
  gnp1max=max(gnp1);
  if gnp1max > 1.01 then
   x_message(['Attenzione   g>1' ; 'Probabile schema di integrazione instabile ']);
   break;
  end;
  xset("window",0)
  xset("wwpc");
  xtape("clear",0);
  xtitle('Evoluzione della soluzione ', 'x', 'g');
  plot2d(X,gnp1,-1,rect=[0.,0.,1.,1.]);
  xset("wshow");
  xset("window",1);
  res=sum(abs((gnp1-gn))/dt)/N;
  plot2d(nn,res,-1,logflag="nl", rect=[0,Resmin,itermax,10.]);
  t=t+dt;
  gn=gnp1;
  nn=nn+1;
end,


xset("window",0);
ge=(exp(lambda*X/mu)-1.0)/(exp(lambda/mu)-1.0);
plot2d(X,ge,2,rect=[0.,0.,1.,1.]);
xset("wshow");


xset("window",1);
//xtape("replay",1); 

//errore numerico
//xset("window",2);
//plot2d(X,ge-gnp1);

//errore di troncamento
//et=ge-L*ge;
//plot2d(X,et);

cont1=x_dialog('Lettura soluzione in un punto?','Y');
if cont1=="Y" then
 cont2="Y";
 while cont2=="Y",
 inode=evstr(x_dialog('Scegli il nodo da valutare(1-N)','0'));
 testo1(1) =strcat( ['Confronto tra la soluzione numerica e la soluzione esatta in']);
 testo1(2) =strcat( ['                    t: ' string(msprintf('%f', t           ))]);
 testo1(3) =strcat( ['                    X: ' string(msprintf('%f', X(inode)    ))]);
 testo1(4) =strcat( ['----------------------------------------------------------------']);
 testo1(5) =strcat( ['soluzione -->  Esatta:' string(msprintf('%f',ge(inode)     )) ' (a stazionario)' ]);
 testo1(6) =strcat( ['             Numerica:' string(msprintf('%f',gnp1(inode)   )) ' (Res:' string(msprintf('%e', res )) ')'  ]);
 testo1(7) =strcat( ['----------------------------------------------------------------']);
 testo1(8) =strcat( ['Soluzione numerica ottenuta con:']);
 testo1(9) =strcat( [ string(msprintf('%s',tschema(nschema))) ]);
 testo1(10) =strcat( ['integrando con  dx:'  string(msprintf('%f',dx )) '      dt:' string(msprintf('%f',dt ))]);
 x_message(testo1);
 cont2=x_dialog('Valutazione di un altro nodo ?','Y');
 end,
end,
cont0=x_dialog('Nuova valutazione ?','Y');
end,


xdel(0);
xdel(1);
quit;