//
// Problema da studiare:
//   .
//   x = A x + g   ,   t in [0,tf]
//
// Condizione iniziale:
//
//   x(0) = x0
//
// L'integrazione si arresta quando x(t) ha valore negativo.
//
clear;
//
A = -1;
//
function y = g(t)
    y = sin(t);
endfunction
//
function y = d1g(t)
    y = cos(t);
endfunction
//
function y = d1x(t,x)
    y = A*x + g(t);
endfunction
//
function y = d2x(t,x)
    y = A*d1x(t,x) + d1g(t);
endfunction
//
function y = g_ev(t,x)
    y = x;
endfunction
//
// INTEGRAZIONE NUMERICA *******************************
//
Percorso = "/home/ciampa/Lavoro/Scilab/Lavoro/LMV/ODE/";
exec(Percorso + 'LMV_TS_1_pv_e.sci');
//
x0 = 4;
t0 = 0;
tf = 10;
EL_MAX = 1d-6;
printf("\nErrore locale massimo consentito = %4.3e\n", EL_MAX);
//
procedura = "LMV_TS_1_pv_e";
[T,U,PASSO,StimaEL] = LMV_TS_1_pv_e(x0,t0,tf,d1x,d2x,g_ev,EL_MAX,"zitto");
//
printf("\nFinito integrazione numerica...\n")
//
// Fine integrazione numerica **************************
//
// Soluzione esatta
//
nota = "si";
if nota == "si" then
//
function y = xe(t)
    y = 9/2 * exp(-t) + (sin(t) - cos(t))/2;
endfunction
//
// errore
//
et = xe(T) - U;
ET = [];
for k = 1:length(T), ET(k) = norm(et(k)); end;
end;
//
// Estremo superiore per l'errore totale
//
SUP = [];
SUP(1) = 0;
for k = 2:length(T),
    SUP(k) = EL_MAX + SUP(k-1)*exp(-PASSO(k-1));  // equazione stabile
end;
//
// grafici ***********************************
//
scf(0); clf();
subplot(321);
if nota == "si" then
   plot(T,et,"r");
   legend("et","in_upper_right");
   xtitle("errore totale");
else
   plot(T,U,"r");
end;
xgrid();
//
subplot(322);
plot(T(1:$-1),PASSO,"r");
xgrid();
xtitle("passo");
//
subplot(323);
if nota == "si" then
   plot(T(2:$-2),log10(ET(2:$-2)),"r",T(2:$-2),log10(SUP(2:$-2)),"b",...
        T(2:$-2),log10(StimaEL(2:$-2)),"g");
   legend("ET","SUP","StimaEL");
else
   plot(T(2:$),log10(SUP(2:$)),"b");
   legend("SUP",2);
end;
xtitle("log10 errore");
xgrid();
//
subplot(324);
if nota == "si" then
   tt = T;
   if max(PASSO) > (tf-t0)/100  then
      tt = linspace(t0,T(1,$),100);
   end;
   plot(tt,xe(tt),"r",T(1:$),U(1:$),"b");
   legend("xe","U");
   xtitle("soluzione")
   xgrid();
end;
//
subplot(325);
a = 0; b = 0.2;
xstring(a,b+0.4,"Problema: dx/dt = -Ax + g");
xstring(a,b+0.2,"Procedura: " + procedura);
xstring(a,b+0.1,"SUP(k) = EL_MAX + SUP(k-1)*exp(-PASSO(k-1))");
format("e");
xstring(a,b+0.0,"EL_MAX = " + string(EL_MAX));
//
subplot(326);
if nota == "si" then
   xstring(a,b+0.3,"Errore totale massimo = " + string(max(ET)));
end;
format("v");
  xstring(a,b+0.2,"Numero passi = " + string(length(PASSO)));
format("e");
//
