%% Finite differences for the advection-diffusion equation % % - D u_xx + v u_x = 1 % % on the domain x=[0,1], % with boundary conditions u(0)=u(1)=0. % % Computational Science and Engineering % Chapter 1, Section 1.2, p 25, Question 19. % clear clc close all % try experimenting with parameters D, v, and N. % D=1; v=1; D=0.05; v=0.2; % D=0.05; v=2; % boundary layer N=5+1; % mesh points xL=0; xR=1; h=(xR-xL)/(N-1); % h is the spatial discretization size x=xL:h:xR; x=x'; % True solution: r=v/D; C=1/(v*(exp(r)-1)); uTrue=x/v+C*(1-exp(r*x)); % diffusion matrix e = ones(N,1); K = spdiags([e -2*e e], -1:1, N, N); K=D*K/h^2; %K=full(K) % uncomment this line to see the matrix K (it will be easier if N is small) % centered difference matrix Del0 = spdiags([-e e], [-1 1], N, N); Del0=v*Del0/(2*h); %Del0=full(Del0) % upwind difference matrix DelPlus = spdiags([-e e], 0:1, N, N); DelPlus=v*DelPlus/h; %DelPlus=full(DelPlus) %% Linear solve A=Del0-K; % centered advection matrix, plus diffusion matrix A(1,:)=0; A(end,:)=0; A(1,1)=1; A(end,end)=1; % boundary conditions b=ones(size(x)); b(1)=0; b(end)=0; % boundary conditions u=A\b; % solve with backslash (centered) A=DelPlus-K; % upwind advection matrix, plus diffusion matrix A(1,:)=0; A(end,:)=0; A(1,1)=1; A(end,end)=1; % boundary conditions U=A\b; % solve with backslash (upwind) %% plot results FS='FontSize'; figure(1) plot(x,u, x,U, '--', x, uTrue, '+', 'LineWidth', 2) title({'- D u_{xx} + v u_x = 1 ', [' D=', num2str(D), ' v=', num2str(v), ' h=', num2str(h)]}, FS, 16) legend('u centered', 'U upwind', 'true solution') xlabel('x', FS, 16) ylabel('u', FS, 16) disp('Finished.')