Lezioni del 13 e 20 aprile 2007

Metodi Numerici per l’equazione del trasporto lineare

Vari dati iniziali: Funzione di Heaviside:

File h0.m

:open:

%
% Funzione di Heaviside (rovesciata)
%
%   -----+
%             |
%             +---------------
%
function y=h0(x)
%
%  se x`e' uno scalare si puo' scrivere cosi:
%
%  if (x > 0)
%    y = 0;   
%  else
%    y = 1;
%  end 
%
%  nel caso vettoriale bisogna fare un ciclo sulle componenti di x
%
  y         = zeros(size(x)) ;
  [n,m] = size(x) ;
  for i=1:n
    for j=1:m
      if ( x(i,j) < 0 )
        y(i,j) = 1;
      end 
    end
  end
end

Onda sinusuidale modulata con una campana:

File p0.m

:open:

%
% onda sinusuidale modulata con una campana
%
%  sin(x) / (1+x^2)
%
function y=p0(x)
  y = sin(x) ./ (1+x.^2 );
end

Funzione scalino in \([-1,1]\):

File q0.m

:open:

%
% Funzione scalino in [-1,1]
%
%            +-----+
%            |     |
%  ----------+     +---------------
%
function y=q0(x)
%
%  se x`e' uno scalare si puo' scrivere cosi:
%
%  if (x < -1)
%    y = 0;   
%  else if ( x > 1 )
%    y = 0;
%  else
%    y = 1;
%  end 
%
%  nel caso vettoriale bisogna fare un ciclo sulle componenti di x
%
  y         = zeros(size(x)) ;
  [n,m] = size(x) ;
  for i=1:n
    for j=1:m
      if ( x(i,j) > -1 && x(i,j) < 1 )
        y(i,j) = 1;
      end 
    end
  end
end

Il seguente programma costruisce la soluzione esatta del problema (con dato iniziale p0.m):

File p.m

:open:

%
% Soluzione esatta del problema
%
% q (x,t) - x t q (x,t) = 0  dato iniziale q(x,0) = q0(x)
%  t             x
%
%
%  q(x,t) = q0(x e^(t^2/2) )
%
function y=p(x,t)
  y = p0( x .* exp(t^2/2) ) ;
end

\[q_t(x,t)-xtq(x,t)=0\]

Il seguente programma costruisce la soluzione esatta del problema (con dato iniziale q0.m):

File q.m

:open:

%
% Soluzione esatta del problema
%
% q (x,t) - x t q (x,t) = 0  dato iniziale q(x,0) = q0(x)
%  t             x
%
%
%  q(x,t) = q0(x e^(t^2/2) )
%
function y=q(x,t)
  y = q0( x .* exp(t^2/2) ) ;
end

Programmi <B>driver</B> per verificare varie soluzioni:

File test.m

:open:

hold off ;
x = [-2:0.01:2] ;
plot(x,q0(x),'.-b','LineWidth',1.5) ; hold on ;
plot(x,q(x,1),'+-r') ; hold on ;
plot(x,q(x,2),'o-y') ; hold on ;
plot(x,q(x,3),'.-m') ; hold on ;
plot(x,q(x,4),'.-g') ; hold on ;

legend('q0','t=1','t=2','t=3','t=4');
axis([-2,2,-0.1,1.1]);

File test1.m

:open:

hold off ;
x = [-2:0.01:2] ;
plot(x,p0(x),'.-b','LineWidth',1.5) ; hold on ;
plot(x,p(x,1),'+-r') ; hold on ;
plot(x,p(x,2),'o-y') ; hold on ;
plot(x,p(x,3),'.-m') ; hold on ;
plot(x,p(x,4),'.-g') ; hold on ;

legend('p0','t=1','t=2','t=3','t=4');
axis([-2,2,-0.5,0.5]);

File test2.m

:open:

hold off ;
beta = -0.5 ;
x = [-2:0.01:2] ;
plot(x,q0(x),'.-b','LineWidth',1.5) ; hold on ;
plot(x,esatta(x,1,beta,@q0),'+-r') ; hold on ;
plot(x,esatta(x,2,beta,@q0),'o-y') ; hold on ;
plot(x,esatta(x,3,beta,@q0),'.-m') ; hold on ;
plot(x,esatta(x,4,beta,@q0),'.-g') ; hold on ;

legend('p0','t=1','t=2','t=3','t=4');
axis([-2,2,-0.1,1.1]);

File test3.m

:open:

hold off ;
alpha = 0.5 ;
x = [-5:0.01:5] ;
plot(x,h0(x),'.-b','LineWidth',1.5) ; hold on ;
plot(x,trasporto(x,1,alpha,@h0),'+-r') ; hold on ;
plot(x,trasporto(x,2,alpha,@h0),'o-y') ; hold on ;
plot(x,trasporto(x,3,alpha,@h0),'.-m') ; hold on ;
plot(x,trasporto(x,4,alpha,@h0),'.-g') ; hold on ;

legend('p0','t=1','t=2','t=3','t=4');
axis([-5,5,-0.1,1.1]);

Programma driver che permette di generare la soluzione esatta dato un profilo iniziale:

File trasporto.m

:open:

%
% Soluzione esatta del problema
%
% q (x,t) + alpha q (x,t) = 0  dato iniziale q(x,0) = q0(x)
%  t                           x
%
%
%  q(x,t) = q0(x - alpha * t )
%
function y=trasporto(x,t,alpha,q0)
  y = q0( x - alpha .* t ) ;
end

---