x_message(['Integrazione problema n.3  con';...
           'schema BTCS.';'Effetti di addensamento nodi ';...
           ' '; ...
           'Ver 3.0  (19/3/2009)'])


function [s]=strech(N,P,Q)
 AN=N-1;
 DETA=1./AN;
 TQI=1./tanh(Q);
 for L=1:N
  AL=L-1;
  ETA=AL*DETA;
  DUM=Q*(1-ETA);
  DUM=1.-tanh(DUM)*TQI;
  s(L)=P*ETA+(1.-P)*DUM;
 end,
endfunction

lambda=1.0;
mu=0.1;
N=11;
dt=0.499e-1;
itermax=1000;


cambia=2;
while cambia==2,


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


txt=['N. nodi';'dt    ';'N. max iter';'Res. min'];
sig=x_mdialog('Parametri  integrazione ',txt,[string(N);string(dt);string(itermax);'1e-5']);
N=evstr(sig(1));
dt=evstr(sig(2));
itermax=evstr(sig(3));
Resmin=evstr(sig(4));

//N=101;
dxi=1/(N-1);
c=lambda*dt/dxi;
s=mu*dt/dxi^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']);
nschema=2;


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',dxi  ))']);
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,

cambia=2;
while cambia==2,


txt=['P ([1.01 ..1.99] ']
sig=x_mdialog('Parametro di addensamento ',txt,['1.01']);
P=evstr(sig(1));
//P=1.8;
Q=2.0;
X=strech(N,P,Q);


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



for i=2:N-1
  dum= (X(i+1)-X(i-1))/(2*dxi);
  dxidx(i)=1./dum;
end
dxidx(1)=3*dxidx(2)-3*dxidx(3)+dxidx(4);
dxidx(N)=3*dxidx(N-1)-3*dxidx(N-2)+dxidx(N-3);
for i=2:N-1

  dxidx2(i)=(dxidx(i+1)-dxidx(i-1))/(Xi(i+1)-Xi(i-1))*dxidx(i);
end




iter=2000;


  A=zeros(N,N);
  A(1,1)=1.;
  for i=2:N-1
    d=mu*dt/dxi*dxidx2(i);
    A(i,i+1)=c*dxidx(i)/2-s*dxidx(i)^2-d/2.;
    A(i,i)=1+2*s*dxidx(i)^2;
    A(i,i-1)=-(c*dxidx(i)/2+s*dxidx(i)^2)+d/2;
  end
  A(N,N)=1.;
  L=inv(A);



//xbasc();
subplot(311);

xtitle('Evoluzione della soluzione ', 'x', 'g');
plot2d(X,gn,-1,rect=[0.,0.,1.,1.]);
s1=gce();




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("wshow");

  res=sum(abs((gnp1-gn))/dt)/N;

    s1.children.data=[X gnp1];

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


ge=(exp(lambda*X/mu)-1.0)/(exp(lambda/mu)-1.0);
plot2d(X,ge,2,rect=[0.,0.,1.,1.]);
s2=gce();



subplot(312);

//errore numerico

subplot(313);
//xset("window",2);
xtitle('Errore globlale', 'x', 'e');
plot2d(X,ge-gnp1);

//errore di troncamento
subplot(312);
xtitle('Errore di troncamento', 'x', 'Et');
et=(ge-L*ge)/dt;
plot2d(X,et);



testo1(1)='cambiare addensamento';
cambia= x_message(testo1,['No' 'Yes']);
subplot(311);
delete(s1.children);
delete(s2.children);
end,

quit;