sandbox/easystab/diffusion_eigenmodes.m

    Eigenmodes of the diffusion equation

    This is a pedagogical introduction for the computation of the eigenmodes of a 1D system. We then use these eigenmodes to see the evolution in time of an initial condition.

    clear all; clf
    % parameters
    N=20; % number of gridpoints
    L=pi; % domain length
    mu=1; % diffusion coefficient
    

    Theory

    In this program we consider the 1D diffusion equation :

    \displaystyle \frac{\partial T}{\partial t} = \mu \frac{\partial^2 T}{\partial x^2} \qquad \textrm{ for } \quad x \in [0,L] with boundary conditions : \displaystyle T_{x=0} = T_{x=L} = 0

    We consider solutions in eigenmode form : \displaystyle T = \hat{T} e^{\lambda t}

    Eigenpairs (\lambda, \hat{T}) are solutions of the following problem :

    \displaystyle \lambda \hat{ T } = \mu \frac{\partial^2 \hat{T}}{\partial x^2}

    The analytical solution is as follows : \displaystyle \lambda_n = - \frac{\mu n^2 \pi^2}{L^2} ; \quad \hat{T}_n = \sin ( n \pi x /L )

    stheory = -mu*pi^2/L^2*[1:1:N].^2;
    

    Numerical resolution

    Discretization

    We use the function dif1D from the easystab project to construct the grid x and the differentiation matrices dx and dxx (the ‘weight’ variable wx is not used here)

    [dx,dxx,wx,x]=dif1D('cheb',0,L,N,3); % try either 'fds' (finite differences) or  'cheb' (chebyshev)
    
    figure ; spy(dxx) % to see the structure of the matrix
    

    Construction of the matrices

    Z=zeros(N,N); I=eye(N); 
    B=I;
    A=mu*dxx; 
    
    % boundary conditions
    loc=[1,N];
    B(loc,:)=0; 
    A(loc,:)=-I(loc,:);
    

    Resolution of the eigenvalue problem

    We compute the eigenmodes using the function eig. We then sort the modes according to decaying real part of the eigenvalue. With this choice, the first eigenvalue will be the one with the largest real part. We then remove the eigenmodes for which the eigenvalue is larger than 1000. We do this because since we have the matrix B to impose the boundary conditions, the system is algebro-differential, withthe concequence that there will be some infinite eigenvalues corresponding to the fact that the constraints are imposed infinitely fast (their dynamics is infinitely rapid).

    % computing eigenmodes
    [U,S]=eig(A,B);
    
    % sort the eigenmodes
    s=diag(S);  [t,o]=sort(-real(s)); 
    s=s(o); U=U(:,o);
    
    
    % show the eigenvalues/eigenmodes
    figure(1);hold off;
    subplot(2,1,1);
    plot(real(stheory),imag(stheory),'ro');
    hold on;
    plot(real(s),imag(s),'bx');
    
    
    xlim([-25,1])
    xlabel('real part');ylabel('imaginary part');title('eigenvalues');
    legend('theory','computed')
    grid on; 
    
    % show the eigenvectors
    subplot(2,1,2);
    co='brkmc';
    for ind=1:3
    plot(x,real(U(:,ind)),co(ind),x,imag(U(:,ind)),[co(ind) '--']);
    hold on
    end
    xlabel('x');  ylabel('eigenvectors');title('eigenvectors');
    grid on; 
    
    set(gcf,'paperpositionmode','auto');
    print('-dsvg','-r80','diffusion_eigenmodes1.svg');
    
    Figure : The eigenvalues and eigenvectors

    Figure : The eigenvalues and eigenvectors

    Time evolution

    Here we use the eigenmodes to show the time evolution of an initial condition. If the initial condition of our system is one of its eigenmodes, we know exactly all the time evolution, it will simply be the eigenvector multiplied by the exponential of the associated eigenvalue time the time. Going further one step, this means that if the initial condition is a combination of the eigenvectors, we get the evolution as sum sum of the eigenvectors weighted by their individual exponential time dependency.

    With the initial condition \displaystyle T(x,0)=\sum_i \alpha_i \hat{T}_i(x) with \alpha_i the amplitude of each eigenmode in the initial condition. The evolution in time is thus \displaystyle T(x,t)=\sum_i \alpha_i \exp(\lambda_i t) \hat{T}_i(x)

    % show the evolution of an initial condition
    figure(2)
    n=5; % number of eigenmodes in the initial condition
    

    Here we simply build a random initial condition by combining the n least stable eigenvectors

    a=randn(n,1); % the weights
    
    
    % time loop
    for t=linspace(0,2,101);
        
        % the present state
        q=U(:,1:n)*(a.*exp(s(1:n)*t));
        plot(x,q);
        hold on
        
        % draw each of the components
        for gre=1:n
           plot(x,U(:,gre)*a(gre)*exp(s(gre)*t),'r--');
        end
        hold off
        axis([0,L,-2,2]); grid on
        title(['t = ',num2str(t)])
        drawnow;
        pause(0.1)
        if (t==0.5)
              set(gcf,'paperpositionmode','auto');
              print('-dsvg','-r80','diffusion_eigenmodes2.svg');
        end      
    end
    

    Here we see in blue the state at the time of the snapshot, and in red the eigenvectors that we have used when building the initial condition, each multiplied with exp(s_i t) where s_i is the associated eigenvalue. This shows that quickly, the evolution of the state of the system tends to the evolution of the least stable eigenmode.

    A snapshot of the time evolution

    A snapshot of the time evolution

    Exercises/contributions

    • Please do the same for other model 1D equations (advection, advection-diffusion…)
    • Please do this for an unstable system