Solving the diffusion equation finite difference technique

[robby991]

Hi everyone. I developed a Matlab program to solve the diffusion equation, a partial differential equation, using the finite difference method. I solved it first with boundary conditions of C = 0, at x=0 and C = 0, at x = 1. I now want to change the boundary condition C = Co (some constant) at x =1. The plot doesn't look right, I need help with the theory because I don't think I quite understand boundry value problems. The Matlab Code is fine. My Matlab code is as follows, and the resulting plot is attached:

numx = 101;   %number of grid points in x
numt = 2000;  %number of time steps to be iterated over 
dx = 1/(numx - 1);
dt = 0.00005;
Co = 3;
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
C(1,1) = 0;    %C=0 at x=0  1st row, 1st column
C(1,numx) = Co; %C=0 at x=1  1st row, numx column
mu = 0.5;
sigma = 0.05;
for i=2:numx-1
   C(i,1) = exp(-(x(i)-mu)^2/(2*sigma^2)) / sqrt(2*pi*sigma^2);
end

%iterate difference equation - note that C(1,j) and C(numx,j) always remain 0
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
end

figure(1);
hold on;
plot(x,C(:,101));
xlabel('x');
ylabel('c(x,t)');

 

[gtoguy97]

title('Solution of Diffusion Equation with C(0)= 0, C(1) = 3');</code>Any help would be greatly appreciated. Thanks!A:The boundary condition C=Co at x=1 suggests that the diffusion equation is solved using Neumann boundary conditions, which are stated by the partial derivative of C with respect to x at x=1.In this case, the finite difference approximation of the partial derivative should be set equal to Co.The formula for a centered difference approximation of the first derivative is<code>dC/dx ~ (C(i+1)-C(i-1))/(2*dx).</code>For the Neumann boundary condition, this should be set equal to Co.<code>(C(numx)-C(numx-1))/(2*dx) = Co.</code>Solving this equation for C(numx) yields<code>C(numx) = Co*2*dx + C(numx-1).</code>This expression can be used to update your code, replacing the boundary condition<code>C(1,numx) = Co;</code>with <code>C(1,numx) = Co*2*dx + C(numx-1);</code>