function ex
%
% Approx mediante elementi finiti del primo ordine
%            - u'' = f, x in (0,2H),    
%            u(0)=Tinf,  u(2H) = Tsup
%
% esercitazione del 12/03/2015

% Dati del problema
Tinf=273; Tsup=293;
kappa=0.6; mu=0.14; H = 1; V=10;
x0=0; x1=2*H;

a=4*V^2*mu/(H^4*kappa);
f=@(x)(a*(H-x).^2);
x=linspace(x0,x1,1000);
u_ex=@(x)(-a/12*(H-x).^4 + (Tsup-Tinf)/(2*H)*x + Tinf + a*H^2/12);   % soluz in forma chiusa

figure(1)
subplot(1,2,1),
plot(x,u_ex(x),'r')
title(' f = f_I')
hold on
figure(2)
subplot(1,2,1),
plot(x,u_ex(x),'r')
title(' (f,v) = quadrature ')
hold on
err_I=[];
err_q=[];
ii=0;

fprintf('\n')
fprintf('   n.raff.    a(e,e) con f_I   a(e,e) con quadr\n');
% loop sul raffinamento della griglia
for k=4:12,

   N = 2^k;
   h = (x1-x0)/N;
   ii = ii+1;
   x = linspace(x0,x1,N+1)'; x_in=x(2:N);

% Crea le nuove matrici di stiffness e di massa
   K = 1/h*spdiags([-ones(N-1,1),2*ones(N-1,1),-ones(N-1,1)],-1:1,N-1,N-1);
   M = h/6*spdiags([ones(N-1,1),4*ones(N-1,1),ones(N-1,1)],-1:1,N-1,N-1);

% Crea il termine noto mediante interpolazione
   F_I = M*f(x_in); 

% Aggiungi al termine noto le condizioni di Dirichlet
   F_I(1) = F_I(1)+Tinf/h;
   F_I(N-1) = F_I(N-1)+Tsup/h;

% Determina i coeff di u_h ed includi i valori al bordo
   alpha = K\F_I;
   uh=[Tinf;alpha;Tsup];

   err_I=[err_I, (u_ex(x_in)-alpha)'*K*(u_ex(x_in)-alpha)];   % norma energia

% Aggiungi al grafico
   figure(1)
   subplot(1,2,1)
   plot(x,uh,'--')
   subplot(1,2,2)
   plot(x,u_ex(x)-uh,'--')

% Crea il termine noto mediante formula di quadratura
   tol=h^2;     % tolleranza richiesta per la quadratura
   Fq=zeros(N-1,1);
   j=1;
   Fq(j) = quad( @(x)(f(x).*(x/h+1-j)),0,h,tol);
   for j=2:N-1
      Fq(j-1) = Fq(j-1)+ quad( @(x)(f(x).*(-x/h+j)),(j-1)*h,j*h,tol);
      Fq(j) = Fq(j)+ quad( @(x)(f(x).*(x/h+1-j)),(j-1)*h,j*h,tol);
   end
   j=N;
   Fq(j-1) = Fq(j-1)+ quad( @(x)(f(x).*(-x/h+j)),(j-1)*h,j*h,tol);

% Aggiungi al termine noto le condizioni di Dirichlet
   Fq(1) = Fq(1)+Tinf/h;
   Fq(N-1) = Fq(N-1)+Tsup/h;
   alphaq = K\Fq;
   uhq=[Tinf;alphaq;Tsup];
   err_q=[err_q, (u_ex(x_in)-alphaq)'*K*(u_ex(x_in)-alphaq)];  %norma energia

   figure(2)
   subplot(1,2,1)
   plot(x,uh,'--')
   subplot(1,2,2)
   plot(x,u_ex(x)-uhq,'--')
   disp(full([k, err_I(ii), err_q(ii)]))
end

fprintf('\n')
fprintf('soluzione esatta nei nodi nel caso di formula di quadratura!\n')

%Ordine di convergenza
fprintf('\n')
fprintf('ordine di convergenza con f_I\n') 
disp(full([log2( err_I(1:end-1)./err_I(2:end))']))

%% risultati
%%>> ex_120315
%%
%%   n.raff.    a(e,e) con f_I   a(e,e) con quadr
%%   4.0000e+00   7.6291e-01   2.0680e-25
%%   5.0000e+00   1.0630e-01   2.9985e-24
%%   6.0000e+00   1.4014e-02   2.7090e-23
%%   7.0000e+00   1.7985e-03   8.5200e-23
%%   8.0000e+00   2.2777e-04   5.2361e-22
%%   9.0000e+00   2.8658e-05   2.5179e-21
%%   1.0000e+01   3.5940e-06   3.1393e-20
%%   1.1000e+01   4.4998e-07   1.5906e-19
%%   1.2000e+01   5.6293e-08   4.7985e-18
%%
%%soluzione esatta nei nodi nel caso di formula di quadratura!
%%
%%ordine di convergenza con f_I
%%   2.8433e+00
%%   2.9232e+00
%%   2.9620e+00
%%   2.9811e+00
%%   2.9906e+00
%%   2.9953e+00
%%   2.9977e+00
%%   2.9988e+00
