Let's consider the diffusion equation with no-flux boundary conditions, and
the initial condition
It is readily seen that this initial condition satisfies the no-flux boundary
conditions, and also that
The following Matlab code solves the diffusion equation for no-flux boundary conditions and for this initial condition:
numx = 101; %number of grid points in x numt = 5000; %number of time steps to be iterated dx = 1/(numx - 1); dt = 0.00005; x = 0:dx:1; %vector of x values, to be used for plotting C = zeros(numx,numt); %initialize everything to zero %specify initial conditions t(1) = 0; %t=0 mu = 0.5; sigma = 0.05; for i=1:numx C(i,1) = 1 - 0.2*cos(pi*x(i)) + 0.1*cos(2*pi*x(i)) + 0.3*cos(3*pi*x(i)) + 0.1*cos(8*pi*x(i)) + 0.05*cos(21*pi*x(i)); end %iterate difference equations for j=1:numt t(j+1) = t(j) + dt; for i=2:numx-1 C(i,j+1) = C(i,j) + (dt/dx^2)*(C(i+1,j) - 2*C(i,j) + C(i-1,j)); end C(1,j+1) = C(2,j+1); %C(1,j+1) found from no-flux condition C(numx,j+1) = C(numx-1,j+1); %C(numx,j+1) found from no-flux condition end figure(1); hold on; subplot(2,3,1); plot(x,C(:,1)); axis([0 1 0.4 1.6]); xlabel('x'); ylabel('c(x,t)'); subplot(2,3,2); plot(x,C(:,11)); axis([0 1 0.4 1.6]); xlabel('x'); ylabel('c(x,t)'); subplot(2,3,3); plot(x,C(:,101)); axis([0 1 0.4 1.6]); xlabel('x'); ylabel('c(x,t)'); subplot(2,3,4); plot(x,C(:,1001)); axis([0 1 0.8 1.2]); xlabel('x'); ylabel('c(x,t)'); subplot(2,3,5); plot(x,C(:,2001)); axis([0 1 0.8 1.2]); xlabel('x'); ylabel('c(x,t)'); subplot(2,3,6); plot(x,C(:,5001)); axis([0 1 0.8 1.2]); xlabel('x'); ylabel('c(x,t)'); %calculate approximation to the integral of c from x=0 to x=1 for j=1:numt+1 s(j) = sum(C(1:numx-1,j))*dx; end figure(2); plot(t,s); xlabel('t'); ylabel('c_{total}'); axis([0 0.25 0.9 1.1]);
We see that the high wavenumber (short wavelength) components of get damped out very quickly - this is what we mean when we say that diffusion smooths out a solution. Notice that by , the only components from the initial condition that are apparent are the ``constant'' component and the ``'' component. The higher the wavenumber (i.e., the shorter the wavelength), the faster the mode is damped out.