//
// funzioni di servizio
//

function [rete,h] = mesh1D (a,b,N)
//costruzione della mesh
h = (b-a)/N;
rete = zeros(1,N+1);
for i=1:N+1
    rete(i) = a+(i-1)*h;
end
endfunction


function [Ki,Hi,Mi] = matrices1D (h)
// Definizione delle matrici locali
Ki = [1 -1; -1 1]/h;
Hi= [-1 1; -1 1]/2;
Mi = h*[2 1; 1 2]/6;
endfunction


function [K,H,M] = assemble1D (Ki,Hi,Mi,N)
// Definizione delle matrici globali

 K = zeros(N+1,N+1);
 H = zeros(N+1,N+1);
 M = zeros(N+1,N+1);

for i=1:N
    K(i:i+1,i:i+1) = K(i:i+1,i:i+1) + Ki;
    H(i:i+1,i:i+1) = H(i:i+1,i:i+1) + Hi;
    M(i:i+1,i:i+1) = M(i:i+1,i:i+1) + Mi;
end
endfunction


function [fi] = localf (rete,h,beta,gamma,u_a,u_b,a,b)
//
// vettore con i valori di f 
//(opportunamente corretta se non omogeneo) sulla mesh
//
  retta=u_a+(rete-a)*(u_b-u_a)/(b-a);
  der_retta=rete*0+(u_b-u_a)/(b-a);
  fi = (h/2)*(f(rete)-beta*der_retta-gamma*retta)';
endfunction


function [rhs] = assemblef (fi,N)
// Costruzione del termine noto f

rhs=zeros(N+1,1);
for i=1:N
    rhs(i:i+1,1) = rhs(i:i+1,1) + fi(i:i+1,1);
end
endfunction

//---------------------------------------------------------
function [value] = uhP1(t, uh, n,a,b)
// valutazione della lineare a tratti  u_h sui punti t
h=(b-a)/n;
for i=1:length(t)
    k=floor((t(i)-a)/h)+1;
    if abs(t(i)-b)<.0000001
        k=n;
    end
    uelem(1) = uh(k);
    uelem(2) = uh(k+1);
    t_loc=(t(i)-(a+(k-1)*h))/h;
    b1=1-t_loc;
    b2=t_loc;
    value(i) = (uelem)' * ([b1, b2])';
end
endfunction
//--------------------------------------------------------------
function [err] = Errore_L2 (a,b,h,uh,N)
// calcolo dell'errore in norma L2  (con trapezi) 
h_tab=h/2;
tab_err=[a:h_tab:b];
uesa=usol(tab_err);
ucalc=uhP1(tab_err,uh,N,a,b);
err=0;
for i=1:max(size(tab_err))-1
    err= err+((uesa(i)-ucalc(i))^2+(uesa(i+1)-ucalc(i+1))^2)*h_tab/2;
end
err=err^(.5);
endfunction
//--------------------------------------------------------------------------------
// FUNZIONI da CAMBIARE a SECONDA del PROBLEMA
//--------------------------------------------------------------------------------
function [value] = f(t)
// valutazione del termine noto f
  value=-4*%pi*cos(2*%pi*t)+(4*(%pi^2)+1)*t.*sin(2*%pi*t);
  //value=t*0;
endfunction
//----------------------------------------------------------------
function [value] = usol(t)
// valutazione della corrispondente soluzione esatta
//  del problema con  parametri alfa=1, beta=0, gamma=1, 
 value = t.*sin(2*%pi*t);
 
 //problema  con parametri: alpha da assegnare,beta=1, gamma=0; 
 //alpha=1/100;
 //app=1/(exp(1/alpha) -1);
 //value=(app*(-1+exp(t/alpha)));
  endfunction 
//---------------------------------------------------------------    
//Soluzione problema ellittico omogeneo unidimensionale
//----------------------------------------------------------------
// -alfa u''+beta u' +gamma u=f in (a,b)
// u(a)=u_a(b)=u_b;
//----------------------------------------------------------------
//
// si immagina di conoscere la soluzione esatta
//
//-----------------------------------------------------------------------
function [uh,h_err,rete,errore,retta] = elliptic1D (a,b,u_a,u_b,N,alfa,beta,gamma)

//dati problema 
//a=0; b=1; u_a=0; u_b=1;
//alfa=1; beta=0; gamma=1;

// discretizzazione
//N=input("assegnare numero di sottointervalli: ");

// costruzione della mesh
[rete,h] = mesh1D(a,b,N);
//
// costruzione della matrice del sistema linerae 
[Ki,Hi,Mi] = matrices1D(h);
[K,H,M] = assemble1D(Ki,Hi,Mi,N);
Atilde = alfa*K +beta*H+gamma*M;
fine=max(size(Atilde));
A = Atilde(2:fine-1,2:fine-1);
//num_cond=cond(A);
//
// costruzione del termine noto f
fi = localf(rete,h,beta,gamma,u_a,u_b,a,b);
rhs = assemblef (fi,N);
rhs = rhs(2:fine-1,1);

//  risoluzione del sistema lineare
uh = A\rhs;
// imposizione delle condizioni al bordo  Dirichlet 

uh = [0; uh; 0];
retta=(u_a+(rete-a)*(u_b-u_a)/(b-a))';
uh=uh+retta;
//
// Grafico della soluzione calcolata e di quella esatta
//
clf(1)
tt=[a:h/10:b];
figure(1)
title('plot of the exact (red) and approximate (blue) solution')
plot(tt,uhP1(tt,uh,N,a,b),tt,usol(tt), 'r')
//
clf(2)
figure(2)
tt=[a:h/3:b];
title('plot of different values of the approximate solution')
plot(tt,uhP1(tt,uh,N,a,b),'*c',rete,uh,'b')
//
// calcolo dell'errore in norma L2 (con trapezi) se la soluzione esatta e' nota
errore=Errore_L2(a,b,h,uh,N);
h_err=h;
endfunction