% This Matlab script calculates the eigenvalues of % the one-dimensional convection equation discretized by % finite differences. The discretization uses central % differences in space and forward Euler in time. % % Periodic bcs are set if periodic_flag == 1. % % Otherwise, an inflow (dirichlet) bc is set and at % the outflow a one-sided (backwards) difference is used. % clear all; periodic_flag = 1; % Set-up grid xL = -4; xR = 4; Nx = 11; % 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 timestep CFL = 1; dt = CFL*dx/abs(u); % Set bc state at left (assumes u>0) UL = exp(-xL^2); % 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); A(i,i+1) = -u/(2*dx); end if (periodic_flag == 1), % Periodic bcs A(1 ,2 ) = -u/(2*dx); A(1 ,Nx-1) = u/(2*dx); A(Nx-1,1 ) = -u/(2*dx); A(Nx-1,Nx-2) = u/(2*dx); else % non-periodic bc's % At the first interior node, the i-1 value is known (UL). % So, only the i+1 location needs to be set in A. A(1,2) = -u/(2*dx); % Outflow boundary uses backward difference A(Nx-1,Nx-2) = u/dx; A(Nx-1,Nx-1) = -u/dx; end % Calculate eigenvalues of A lambda = eig(A); % Plot lambda*dt plot(lambda*dt,'*'); 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;