% Advection diffusion with small viscosity % %%%%%%%%%%%%%%%%%%%%%%%%% % % visc* u_xx - u_x = f % % u(0)=alpha, u(1)=beta, alpha=1, beta=3 % % f=-1 % % In this case (you can verify that) the exact solution is % % u=alpha+x+(beta-alpha-1)*(exp(x/visc)-1)./(exp(1/visc)-1); % %%%%%%%%%%%%%%%%%%%%%%%% % % The solution of the problem is plotted against the numerical % approximation (see below). % %%%%%%%%%%%%%%%%%%%%%%%%% % % You can run the program with different values of the viscosity parameter % visc = 0.3, 0.1, 0.05 % and see that the "boundary layer" gets steaper and steaper. This fact % creates difficulties for the numerical method. % %%%%%%%%%%%%%%%%%%%%%%%%%% % Try also to see what happens if you replace the backward difference % operator D with central differences and modify the boundary term so that % you get an overall second order method in h (the space discretization % parameter). You'll see a bad oscillation near the boundary layer and only % using h small enough (m bigger) this undesired effect will go away. % % This is related to what happens to u_t+au_x=0 (a=1): using central differences % to approximate u_x and forward differences to approximate u_t you get a % method that is always unstable (independently on h and k). % But to fix it you can add a term corresponding to the central difference % discretization of u_xx ("add numerical viscosity") and get the Lax-Wendroff % method (see ch. 7.3). % % You can also try to solve % u_t+au_x=visc* u_xx-f % a=1, u(0)=alpha, u(1)=beta, u(x,0)=g(x) % Using similar methods. You'll get a solution evolving from the initial % function g(x) to the solution of visc*u_xx-au_x=f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % %%%%%%%%%%%%%%%%%%%%%%%%%%% % T=1; alpha=1; beta=3; % m=100; h=T/(m+1); xx=[h:h:T-h]'; visc=0.005; u=alpha+xx+(beta-alpha-1)*(exp(xx/visc)-1)./(exp(1/visc)-1); plot([0;xx;T],[alpha;u;beta]); e=ones(m,1); A=1/(h^2)*spdiags([e -2*e e],-1:1,m,m); %D=1/(2*h)*spdiags([-e e],[-1 1],m,m); D=1/h*spdiags([-e e],-1:0,m,m); F=-ones(m,1); F(1,1)=F(1,1)-alpha*visc/(h^2)- alpha/(h);%alpha/(2*h); F(end,1)=F(end,1)-beta*visc/(h^2); Am=visc*A-D; U=Am\F; hold on plot([0;xx;T],[alpha;U;beta],'r') hold off