function mit18086_fd_transport_growth(method) %MIT18086_FD_TRANSPORT_GROWTH % Applies various finite difference methods for the transport % equation u_t=u_x. Additionally shows the von Neumann growth % factor in complex plane. % Choose method by giving number as parameter. Methods of % choice are (1) upwind, (2) downwind, (3) centered, % (4) Lax-Friedrichs, (5) Lax-Wendroff, and (6) Crank-Nicolson. % 03/2007 by Benjamin Seibold % http://www-math.mit.edu/~seibold/ % Feel free to modify for teaching and learning. n = 90; % number of space gridpoints dt = 1e-2; % time step tf = 2.4; % final time if nargin<1, method = 1; end x = linspace(-1,1,n)'; h = x(2)-x(1); r = dt/h; disp(sprintf('Courant number: %0.2f',r)) u = ic(x); u0 = u; I = eye(n); R = diag(ones(1,n-1),1); Dxc = (R+diag([-1 zeros(1,n-2) 1])-R')/(2*h); Dxx = (R-diag([1 2*ones(1,n-2) 1])+R')/h^2; switch method case 1 name = 'upwind'; A = (R-I)/h; case 2 name = 'downwind'; A = (I-R')/h; case 3 name = 'centered'; A = Dxc; case 4 name = 'Lax-Friedrichs'; A = Dxc+(h^2/dt/2)*Dxx; case 5 name = 'Lax-Wendroff'; A = Dxc+(dt/2)*Dxx; case 6 name = 'Crank-Nicolson'; M = (I-(dt/2)*Dxc)\(I+(dt/2)*Dxc); end if not(exist('M')), M = I+dt*A; end t = linspace(0,2*pi,100); nh = floor(size(M,1)/2)-1; g = M(nh+1,1:2*nh+1)*exp(i*(-nh:nh)'*t); for tn = 1:ceil(tf/dt) u = M*u; clf subplot(1,3,1:2) plot(x,u0,'b:',x,u,'r.-') axis([-1 1 -.1 1.1]) title(sprintf('%s , t=%0.2f',name,tn*dt)) subplot(1,3,3) patch(cos(t),sin(t),.9*[1 1 1]) hold on, plot(g,'r-','linewidth',2), hold off axis equal axis(1.5*[-1 1 -1 1]) title('growth factor') drawnow end function y = ic(x) % initial condition function y = max((1-4*(x-.5).^2),0).^2;