% This Matlab script calculates the eigenvalues of the one-dimensional % convection-diffusion equation discretized by finite differences. % The discretization uses central differences in space and forward % Euler in time. % % Periodic bcs are set. % clear all; % Set-up grid xL = -4; xR = 4; Nx = 21; % number of points x = linspace(xL,xR,Nx); % Calculate cell size in control volumes (assumed equal) dx = x(2) - x(1); % Set velocity u = 1; % Set viscosity using Reynolds based on dx Re = 1e+2; % Don't set to zero identically... use a small positive value mu = u*dx/Re; % Set timestep CFL = 0.5; dt = CFL*dx/abs(u); % Allocate matrix to hold stiffness matrix (A). % A = zeros(Nx-1,Nx-1); % Construct A except for first and last row for i = 2:Nx-2, A(i,i-1) = u/(2*dx) + mu/dx^2; A(i,i ) = -2*mu/dx^2; A(i,i+1) = -u/(2*dx) + mu/dx^2; end % Periodic bcs A(1 ,Nx-1) = u/(2*dx) + mu/dx^2; A(1 ,1 ) = -2*mu/dx^2; A(1 ,2 ) = -u/(2*dx) + mu/dx^2; A(Nx-1,Nx-2) = u/(2*dx) + mu/dx^2; A(Nx-1,Nx-1) = -2*mu/dx^2; A(Nx-1,1 ) = -u/(2*dx) + mu/dx^2; % Calculate eigenvalues of A lambda = eig(A); % Plot lambda*dt plot(lambda*dt,'r*'); xlabel('Real \lambda\Delta t'); ylabel('Imag \lambda\Delta t'); % Overlay Forward Euler stability region th = linspace(0,2*pi,101); hold on; plot(-1 + sin(th),cos(th)); hold off; axis('equal'); grid on;