btn=x_message(["Approssimazione di derivate con";...
             "formule alle differenze finite.";" ";...
            "Calcolo derivate ed analisi degli errori";" "; ...
            "Ver 3.0  (12/3/2009)"])

xset("window",0);

pi=atan (1)*4;

n=x_choose(['y=sin(x)';'y=x*cos(x^2)'; 'y=-x^2+x'],['Scegli la funzione ']);
select n,
 case 1 then function y=f(x),y=sin(x*pi), endfunction;  
 case 2 then function y=f(x),y=x.*cos((x*pi)^2), endfunction,
 case 3 then function y=f(x),y=-x^2+x, endfunction,
end

select n,
 case 1 then function y=df(x),y=pi*cos(x*pi), endfunction,
 case 2 then function y=df(x),y=cos((x*pi)^2)-2*(pi*x)^2.*sin((x*pi)^2), endfunction,
 case 3 then function y=df(x),y=-2*x+1, endfunction,
end




xx=0:1/1000:1;
xbasc();
subplot(121);
plot2d(xx, f(xx),1,"121");

x=evstr(x_dialog('inserire valore di x [0,1]','0.0'));



cont1="Y";
while cont1=="Y"

dx=0.1;
n=x_choose(['I ord. avanti';'I ord. indietro ';'II centr.';'II ord.  avanti'],['Scegli al formula alle differenze finite']);

cont="Y";
ii=1;
while cont=="Y"

subplot(121);
//xbasc();

plot2d(xx, f(xx),1,"121");
xtitle('funzione');
dx=evstr(x_dialog('inserire valore di dx [0,1]',string(dx)));




//n=x_choose(['I ord. avanti';'I ord. indietro ';'II centr.';'II ord.  avanti'],['Scegli al formula alle differenze finite']);

for i=1:3
select n,
 case 1 then dfn(i)=(f(x+dx)-f(x))/(dx);,
 case 2 then dfn(i)=(f(x)-f(x-dx))/(dx);,
 case 3 then dfn(i)=(f(x+dx)-f(x-dx))/(2*dx);,
 case 4 then dfn(i)=(-3*f(x)+4*f(x+dx)-f(x+2*dx))/(2*dx);,
end 
derr(i)=abs(df(x)-dfn(i));
dxx(i)=dx;
dx=dx/2;
end
dx=dx*2;

select n,
 case 1 then xxx=[x x+dxx(1)];            yyy=f(xxx);,
 case 2 then xxx=[x-dxx(1) x];            yyy=f(xxx);,
 case 3 then xxx=[x-dxx(1) x+dxx(1)];     yyy=f(xxx);,
 case 4 then xxx=[x x+dxx(1) x+2*dxx(1)]; yyy=f(xxx);,
end


  xset("mark",1,4)
  plot2d(xxx,yyy,-1 ,"101");
  plot2d(xx, [f(x)+df(x)*(xx-x)],1,"101");
  plot2d(xx, [f(x)+dfn(1)*(xx-x)],2,"101");
//plot2d(xx, [f(x)+dfn(2)*(xx-x)],3,"101");
//plot2d(xx, [f(x)+dfn(3)*(xx-x)],4,"101");


//xset("window",1);
////xbasc();
//xtitle('Analisi di Convergenza', 'Dx', 'Errore');
//plot2d(dxx,[derr],style=[-ii],logflag="ll",rect=[0.0001,0.0001,1.0,10.]);
//plot2d(dxx,[derr],style=[ii],logflag="ll",rect=[0.0001,0.0001,1.0,10.]);

testo(1)=strcat( ['Risultato numerico:']);
testo(2)=strcat( ['Valore esatto della derivata   :' string(df(x))]);
testo(3)=strcat( ['Valore calcolato della derivata:' string(dfn(1))]);
testo(4)=strcat( ['Errore                         :' string(df(x)-dfn(1))]);
testo(5)=strcat( [' ' ' ']);
x_message(testo);

subplot(122);
//xset("window",1);
//xbasc();
xtitle('Analisi di Convergenza', 'Dx', 'Errore');
plot2d(dxx,[derr],style=[-ii],logflag="ll",rect=[0.0001,0.0001,1.0,10.]);
plot2d(dxx,[derr],style=[ii],logflag="ll",rect=[0.0001,0.0001,1.0,10.]);


//testo1(1)=strcat( ['Analisi della convergenza']);
//                 1234567   12345678
//testo1(2)=strcat( ['     dx' '     Dfn']);
//testo1(3)=strcat( [string(msprintf('%2.5f',dxx(1))) ' '  string(dfn(1))  ]);
//testo1(4)=strcat( [string(msprintf('%2.5f',dxx(2))) ' '  string(dfn(2))  ]);
//testo1(5)=strcat( [string(msprintf('%2.5f',dxx(3))) ' '  string(dfn(3))  ]);
//x_message(testo1)

tmp=(dfn(1)-dfn(2))/(dfn(2)-dfn(3));
 if tmp > 0 then 
   ns=log(tmp)/log(2);
   k=(dfn(1)-dfn(2))/(dxx(1)^ns-dxx(2)^ns);
   dfs=dfn(3)-k*dxx(3)^ns;
   testo2(1)=strcat( ['Analisi della convergenza']);
   testo2(2)=strcat( ['Stima ordine ed errore ']);
   testo2(3)=strcat( ['valore esat.:' string(msprintf('%2.5f',df(x)))]);
   testo2(4)=strcat( ['valore Estr.:' string(msprintf('%2.5f',dfs))]);
   testo2(5)=strcat( ['Stima ordine:' string(msprintf('%2.5f',ns))]);
//   //                 1234567   12345678   12345678
   testo2(6)=strcat( ["   dx" "       Dfn" "    errS" "    errE"]);
   testo2(7)=strcat( [string(msprintf('%2.5f',dxx(1))) ' '  string(msprintf('%2.5f',dfn(1))) ' ' string(msprintf('%2.5f',dfs-dfn(1))) ' ' string(msprintf('%2.5f',df(x)-dfn(1)))]);
   testo2(8)=strcat( [string(msprintf('%2.5f',dxx(2))) ' '  string(msprintf('%2.5f',dfn(2))) ' ' string(msprintf('%2.5f',dfs-dfn(2))) ' ' string(msprintf('%2.5f',df(x)-dfn(2)))]);
   testo2(9)=strcat( [string(msprintf('%2.5f',dxx(3))) ' '  string(msprintf('%2.5f',dfn(3))) ' ' string(msprintf('%2.5f',dfs-dfn(3))) ' ' string(msprintf('%2.5f',df(x)-dfn(3))) ]);
  x_message(testo2);
 else
   x_message('Convergenza non monotona');
 end

//df(x);
//[dxx dfn ];
//[dxx derr];
//[ 1 1 ;
// 0.5 derr(1)/derr(2) ;
// 0.25 derr(2)/derr(3) ];
//[ (derr(1)-derr(2))/(derr(2)-derr(3)) ];
cont=x_dialog('Calcolare con altro valore di dx','Y');
ii=ii+1;
end

cont1=x_dialog('Calcolare con altra formula alle differenze','Y');
end

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