% This Matlab script solves the one-dimensional convection % equation using a finite volume algorithm. The problem is % assumed to be periodic. % clear all; % Set-up grid xL = -4; xR = 4; Nx = 4*120; % number of control volumes x = linspace(xL,xR,Nx+1); % Calculate midpoint values of x in each control volume xmid = 0.5*(x(1:Nx) + x(2:Nx+1)); % Calculate cell size in control volumes (assumed equal) dx = x(2) - x(1); % Set velocity u = 1; % Set final time tfinal = 1; % Set timestep CFL = 1.5; dt = CFL*dx/abs(u); % Set initial condition to U0 = exp(-x^2) % Note: technically, we should average the initial % distribution in each cell but I chose to just set % the value of U in each control volume to the midpoint % value of U0. U0 = exp(-xmid.^2); U = U0; t = 0; % Loop until t > tfinal while (t < tfinal), Ubc = [U(Nx), U, U(1)]; % This enforces the periodic bc % Calculate the flux at each interface F = 0.5* u *( Ubc(2:Nx+2) + Ubc(1:Nx+1)) ... - 0.5*abs(u)*( Ubc(2:Nx+2) - Ubc(1:Nx+1)); % Calculate residual in each cell R = F(2:Nx+1) - F(1:Nx); % Forward Euler step U = U - (dt/dx)*R; % Increment time t = t + dt; % Plot current solution stairs(x,[U, U(Nx)]); axis([xL, xR, -0.5, 1.5]); grid on; drawnow; end % Get solution at x=1 xloc = 1; uloc_exact = 1; for i = 1:Nx, if ((x(i) <= xloc) & (x(i+1) > xloc)), uloc = U(i); break; end end fprintf('U(x=%f,t=%f) = %f Error = %f\n',xloc,t,uloc,uloc-uloc_exact);