%Script prova metodo FDM
%Conduzione calore - non stazionario - 1D
%
%Scritto da Pavan Andrea - 19/11/2021
clear;
clc;

%% dati input
dx = 0.1;      %discretizzazione spaziale
x = 0:dx:10;       %griglia
x = x(:);
k = 0.6;        %conducibilità termica [W/(m*K)]
rho = 1.0;      %densità [kg/m3]
c = 4.2;        %calore specifico [J/(kg*K)]
Q = @(x,t) cos(2*x).*(rho*c*pi*cos(pi*t)+4*k*sin(pi*t));       %sorgenti calore [W/m3]
u0 = @(x) 0;     %condizioni iniziali [K]
uL = @(t) sin(pi*t);        %condizioni contorno sx [K]
uR = @(t) cos(20)*sin(pi*t);        %condizioni contorno dx [K]
t0 = 0;     %istante iniziale [s]
tf = 4.5;       %istante finale [s]
ue = @(x,t) cos(2*x).*sin(pi*t);       %soluzione esatta [K]


%% metodo Eulero Esplicito
dt = 0.02;      %passo temporale
t = t0:dt:tf;       %istanti
t = t(:);
u = zeros(length(x),length(t));     %soluzione approssimata
err = 0*t;      %errore massimo
u(:,1) = u0(x);
for i=2:length(t)
    u(1,i) = uL(t(i));
    u(end,i) = uR(t(i));
    u(2:end-1,i) = u(2:end-1,i-1) + dt*((k/(rho*c))*(u(1:end-2,i-1)-2*u(2:end-1,i-1)+u(3:end,i-1))/(dx^2) + Q(x(2:end-1),t(i))/(rho*c));
    err(i) = max(abs(u(2:end-1,i)-ue(x(2:end-1),t(i))));

    %grafico soluzione
    figure(1);
    plot(x,u(:,i),'b-o');
    hold on;
    plot(x,ue(x,t(i)),'r-');
    title(['Conduzione calore 1D instazionario (t=' num2str(t(i)) 's)']);
    legend('Soluzione numerica','Soluzione esatta');
    xlabel('x');
    ylabel('u');
    ylim([-1 1])
    hold off;
end

%grafico errore
figure(2);
semilogy(t,err)
title('Andamento errore nel tempo');
xlabel('t');
ylabel('err');

