# 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');


# 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.