
% Equazione delle onde

% coefficienti della PDE:
ro=0.3;
v=0.8;
c=0.5;
sigma=0; % consideriamo solo il caso "sigma=0"

% discretizzazione nello spazio
N=18;     % numero di nodi della griglia
% qui serve anche il passo della griglia (viene usato per verificare la condizione di stabilità di Eulero in avanti)
h=1/(N-1);     % passo della griglia

% costruisco la matrice "A" (approssimazione del laplaciano):
K = mat_ellipt_2D(N,sigma);
% matrice di massa M e di smorzamento G
M=ro*eye(size(K,2));
G=v*eye(size(K,2));

%u_vera=ones(size(A,2),1); %voglio che la soluzione esatta sia questa
%b=K*u_vera; 

x=[h:h:1-h]';
y=[h:h:1-h]';
[xx,yy]=meshgrid(x,y);


% discretizzazione nel tempo
figure(2);clf, title('x(T)');

% ottengo l'equazione differenziale del secondo ordine: M * u"(t)+G * u'(t)+K * u(t)=f(t)
% da questo passo ad un'equazione differenziale del primo ordine, raddoppiando la dimensione di u
u0=zeros(size(K,2),1);
d_u0=ones(size(K,2),1);
x0=[d_u0; u0];

A=[-inv(M)*G  -inv(M)*K;  eye(size(K,2))  zeros(size(K,2))];   % matrice del nuovo sistema di equazioni di dimensione doppia

t0=0;
T=5;
dt = 0.05;
tt=[t0:dt:T];
x=x0;
for t=tt(2:end-1)
    % costruisco un termine forzante oscillante
    f=zeros(size(K,2),1);
    k=1;
    for j=1:(N-2)
       for i=1:(N-2)
           f(k)=sin(x(i)*y(j))-cos(t);
           k=k+1;
       end
    end
    b=[inv(M)*f; zeros(size(K,2),1)];    
    v=x(:,end);
    v = (eye(2*size(K,2)) - dt*A) \ (v + dt*b);
    x=[x v];
    % Grafico della soluzione in "t", riportata alle coordinate della griglia
    sol_T=reshape(x(size(K,2)+1:end,end),N-2,N-2);
    mesh(xx,yy,sol_T); axis([0 1 0 1 -0.15 0.1]), drawnow, tic, 
    while toc<0.05, end;
end






