Facing problem to solve convection-diffusion equation

I want to solve simple diffusion equation MATLAB.

 (del_C/del_t)=(D^2)(del^2_C/del_x^2)-U(del_C/del_x)

I.C: C(x,0)=C0; B.C: C(0,t)=C0, C(L,t)=C1

My program:

numr =50;         %number of grid points in x direction
numt = 10000;       %number of time steps 
L = 0.01;           %thickness of a flat plate (m)
dx = R/(numr - 1);  %grid size in x direction
dt = 0.0005;
D = 1.0e-05;        %species diffusivity 
C0 = 1.0;           %species initial concentration
Cinf = 0.1;         %species surface concentration 
x = 0:dx:L;         %vector of x positions

C = zeros(numr,numt);  %initialize the concentration C(x,t) matrix

%specify initial conditions along the radius(C=C0,x=0 to L,t=0)
t(1) = 0;                 
for i=1:numr
  C(i,1)=C0; 
end

% uniform boundary condition at the centre(C=C0,x=0,t=0 to tf)constant mass
% source term
 for j = 1:numt
     C(1,j)=C0; 
 end

% boundary condition at particle surface (C=Cinf,x=L,t>0)
% specis concentration is fixed at the boundary
for j=2:numt
  C(numr,j)=Cinf; 
end

%Stability criteria
%dx<(2*D)/u
%dt<(dx^2)/(2*D)

u=0.5;

s=D*dt/(2*dx^2);
p=u*dt/(4*dx);

%iteration based on forward differencing scheme 
for j=2:numt
   t(j) = t(j-1) + dt;
   for i=2:numr-1
C(i,j+1)=C(i,j)+(s-p)*(C(i+1,j+1)+C(i+1,j))-(2.0*s)*(C(i,j+1)+C(i,j))+(s+p)*(C(i-1,j+1)+C(i-1,j));
end
end

plot(C,x)
xlabel('x direction')
ylabel('C')

Answer 1

[Walter Roberson]

Your line

C(i,j+1)=C(i,j)+(((s-p)*C(i+1,j+1)+(s-p)*C(i+1,j)+((-2*s)*C(i,j+1))+((-2*s)*C(i,j))+((s+p)*C(i-1,j+1))+((s+p)*C(i-1,j));

is missing two closing brackets somewhere

C(i,j+1)=C(i,j)+(((s-p)*C(i+1,j+1)+(s-p)*C(i+1,j)+((-2*s)*C(i,j+1))+((-2*s)*C(i,j))+((s+p)*C(i-1,j+1))+((s+p)*C(i-1,j));
 1     0  1   0 123   2  3       2 3   2  3     2 34    3  4     32 34    3  4   32 34   3  4       32 34   3  4     32

so at the end of the expression you have 2 unclosed brackets.

[Raj001]

Bracket is not a problem. Mistakenly put here.

After running the program, Index exceeds matrix dimensions is displaying.

[Walter Roberson]

Difficult for me to check without knowing the correct code.

[Raj001]

I have corrected the brackets. Also I have modified the expression. This one is correct code.

Answer 2

[Torsten]

On the right-hand side of the equation

C(i,j+1)=C(i,j)+(((s-p)*C(i+1,j+1)+(s-p)*C(i+1,j)...

you refer to elements C(.,j+1) that are not yet known. This is the general problem.

The "index exceeds matrix dimensions" error is caused by this reference to C(i,j+1) for j=numt since you allocated C as C=zeros(numr,numt).

Best wishes

Torsten.

[Raj001]

Thanks. Matrix dimension now is ok.

[Torsten]

But the results will be incorrect. You have to solve a tridiagonal linear system of equations to get the C(.,j+1) from the C(.,j).

Best wishes

Torsten.