messagebox(['Integrazione problema n.3  con' 'schemi alle differenze finite.' ' '  'Ver 3.1  (24/3/2010)'], 'Program3',"modal", "Ok" );


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)';...
         'FTUS (esplicito,  I ord. tempo, I ord. spazio)';...
         'CN   (implicito, II ord. tempo, II ord. spazio)'];
nschema=x_choose_modeless(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( ['                              c+2s:' string(msprintf('%2.5f',c+2*s ))']);
testo(12)=strcat( ['                       N. iter max:' string(msprintf('%2.5f',itermax))']);
testo(13)=strcat( ['                         Res. min.:' string(msprintf('%2.5f',Resmin))']);
testo(14)=strcat( ['            Schema di integrazione:' string(msprintf('%s',tschema(nschema) ))']);
testo(15)=strcat( ['------------------------------------------------------------']);
testo(16)=strcat( ['Vuoi cambiare i parametri?']);

cambia=messagebox(testo,'Program3',"question",['No' 'Yes'], "modal");
 
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);
elseif nschema==3 then
  //FTUS
  A=zeros(N,N);
  A(1,1)=1.;
  for i=2:N-1
    A(i,i+1)=s;
    A(i,i)=1-c-2*s;
    A(i,i-1)=s+c;
  end
  A(N,N)=1.;
  L=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,




xbasc();
subplot(121);
xtitle('Evoluzione della soluzione ', 'x', 'g');
subplot(122);
xtitle('Storia del residuo ', 'iter.', 'Res.');
subplot(121);
plot2d(X,gn,-1,rect=[0.,0.,1.,1.]);
s1=gce();

//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.5  then
   messagebox(['Attenzione   g>1','Probabile schema di integrazione instabile'], 'Program3',"modal", "Ok" );
   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.]);
    subplot(121);
    s1.children.data=[X gnp1];
    xset("wshow");
    subplot(122);
    res=sum(abs((gnp1-gn))/dt)/N;
//    plot2d(nn,res,-1,logflag="nl", rect=[0,Resmin,itermax,10.]);
    plot2d(nn,res,-2,logflag="nl", rect=[0,Resmin,itermax,10.],frameflag=1);

    t=t+dt;
    gn=gnp1;
    nn=nn+1;
end,


subplot(121);
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 ))]);
 messagebox(testo1,'Risultato',"modal", "Ok");
 cont2=x_dialog('Valutazione di un altro nodo ?','Y');
 end,
end,
cont0=x_dialog('Nuova valutazione ?','Y');
end,


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