% soluzione problema ellittico  unidimensionale
%-------------------------------------------------------------------------
% -alpha 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,mesh,err] = elliptic1D (a,b,u_a,u_b,N,alfa,beta,gamma)

close all
%
% costruzione della mesh
[mesh,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;
A = Atilde(2:end-1,2:end-1);
spy(A)
%
% costruzione del termine noto f
fi = localf(mesh,h,beta,gamma,u_a,u_b,a,b);
rhs = assemblef (fi,N);
rhs = rhs(2:end-1,1);

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

uh = [0; uh; 0];
retta=(u_a+(mesh-a)*(u_b-u_a)/(b-a))';
uh=uh+retta;
% Grafico della soluzione calcolata e di quella esatta

tt=[a:h/10:b];
figure
hold on
title('plot of the exact and approximate solution');
plot(tt,uhP1(tt,uh,N,a,b))
plot(tt,usol(tt), 'r');
hold off

% calcolo dell'errore in norma L2 discreta (con trapezi) 
%se la soluzione esatta e' nota
err=Errore_L2(a,b,h,uh,N);
err=0;
return
%--------------------------------------------------------------------------
% FUNZIONI di SERVIZIO
%--------------------------------------------------------------------------
% creazione della mesh
function [mesh,h] = mesh1D (a,b,N)
h = (b-a)/N;
mesh =zeros(1,N+1);
for i=1:N+1
    mesh(i) = a+(i-1)*h;
end
return

%--------------------------------------------------------------------------
% definizione delle matrici locali
function [Ki,Hi,Mi] = matrices1D (h)
Ki = [1 -1; -1 1]/h;
Hi= [-1 1; -1 1]/2;
Mi = h*[2 1; 1 2]/6;
return
%--------------------------------------------------------------------------
% definizione delle matrici globali
function [K,H,M] = assemble1D (Ki,Hi,Mi,N)
K = sparse(N+1,N+1);
H = sparse(N+1,N+1);
M = sparse(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
return
%--------------------------------------------------------------------------
% valori di f sulla mesh
% con opportuna correzione per problema non omogeneo
function [fi] = localf (mesh,h,beta,gamma,u_a,u_b,a,b)
retta=u_a+(mesh-a)*(u_b-u_a)/(b-a);
der_retta=mesh*0+(u_b-u_a)/(b-a);
fi = (h/2)*(f(mesh)-beta*der_retta-gamma*retta)';
return
%--------------------------------------------------------------------------
% Costruzione del termine noto 
function [rhs] = assemblef (fi,N)
rhs = sparse(N+1,1);
for i=1:N
    rhs(i:i+1,1) = rhs(i:i+1,1) + fi(i:i+1,1);
end
return

%--------------------------------------------------------------------------
% Calcolo dell'errore in norma L2
function [err] = Errore_L2 (a,b,h,uh,N)
% calcolo dell'errore in norma L2 discreta (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);
return
%--------------------------------------------------------------------------
% valutazione  della lineare a tratti u_h in t
function value = uhP1(t, uh, n,a,b)
    h=(b-a)/n;
for i=1:length(t)
    k=floor((t(i)-a)/h)+1;
    if t(i)==b
        k=n;
    end
    uelem(1) = uh(k);
    uelem(2) = uh(k+1);
    t_loc=(t(i)-(a+(k-1)*h))/h;
    [b1, b2]= basisP1(t_loc);
    value(i) = (uelem) * ([b1, b2])';
end
return
%--------------------------------------------------------------------------
function [b1, b2, bx1, bx2]= basisP1(t)
% valutazione locale della base dei polinomi di grado 1
b1 = (1-t);
b2 = t;
% valutazione delle derivate
bx1=-1;
bx2=1;
return
%--------------------------------------------------------------------------
% FUNZIONI da CAMBIARE a seconda del problema
%-------------------------------------------------------------------------
% valutazione del termine noto  f
function value = f(t)
 %esempio con alpha=1, beta=0, gamma=1
 value=-4*pi*cos(2*pi*t)+(4*(pi^2)+1)*t.*sin(2*pi*t);   
return
%--------------------------------------------------------------------------
% valutazione della soluzione esatta
function value = usol(t)    
    value= t.*sin(2*pi*t); 
return