Coupled diffusion equations code

[robby991]

Hi, I am trying to solve a problem using the central difference scheme in Matlab where there are two coupled diffusion equations (where both diffusion equations depend on S(i,j)). They are as follows:

dS/dt = Ds * d^2S/dx^2 - Vmax[B]S[/B]/Km+[B]S[/B]
and 

dP/dt = Dp * d^2P/dx^2 - Vmax[B]S[/B]/Km+[B]S[/B]

with boundary conditions

dS/dx(x=0) = 0;  S(x=d,t) =0;

and

P(x=0,t)=0; P(x=d,t)=0;

I am unsure how to setup the code and I was wondering if anyone can help. From what I understand I need to setup only 1 grid, call it Z, with grid spacing in the space (numt) and time (numx) domain. So first inititalize it to zeros.

Z = zeros(numt,numx);

Now I am confused as to how the loop will be setup for iterations, as there will be two different equations with their associated boundary conditions stepping through the same grid. I don't understand how this can be done without overwriting each other. Here is some code I wrote, where I assumed 2 separate grids, 1 for S and 1 for P. I would appreciate some help in re-writing it for 1 grid, or if someone can provide an explanation on how to do it. Thanks.

%specify initial conditions%

t(1) = 0;                 %1st t position = 0
S(1,:) = S0;                        
P(:,numx)=0;

% Dirchelet Boundary Conditions, where Neumann BC dS/dx=0 is in body of for loop
S(:,numx) = S0;          %S(x=d,t)=S0
P(:,1) = 0;              %P(x=0,t)=0
P(:,numx) = 0;           %P(x=d,t)=0

%iterate central difference equation% 

for j=1:numt-1  
    
      
    %2nd Derivative Central Difference Substrate Iteration% 
     
    for i=2:numx-1
      S(j+1,i) = S(j,i) + (dts/dxs^2)*Ds*(S(j,i+1) - 2*S(j,i) + S(j,i-1))-((Vmax*dts*S(j,i))/(Km+S(j,i))); 
    end
    
    S(j+1,1)=S(j,1)+dts*Ds*2*(S(j,2)-S(j,1))./dxs.^2-((Vmax*dts*S(j,i))/(Km+S(j,i)));      %Neumann Boundary Condition
    
end

for j=1:numt-1  
    
      
    %2nd Derivative Central Difference Iteration% 
     
    for i=2:numx-1
      P(j+1,i) = P(j,i) + (dtp/dxp^2)*Dp*(P(j,i+1) - 2*P(j,i) + P(j,i-1))+((Vmax*dtp*S(j,i))/(Km+S(j,i))); 
    end
    
end

 

[gtoguy97]

The answer is that you can use the same grid for both equations. You should set up a loop that iterates through each of the equations in turn, updating the grid at each step. The boundary conditions should be applied at each step. For example: %specify initial conditions%t(1) = 0; %1st t position = 0Z(:,numx)=0;% Dirchelet Boundary Conditions, where Neumann BC dS/dx=0 is in body of for loopZ(:,1) = 0; %P(x=0,t)=0Z(:,numx) = 0; %P(x=d,t)=0%iterate central difference equation% for j=1:numt-1 %2nd Derivative Central Difference Substrate Iteration% for i=2:numx-1 Z(j+1,i) = Z(j,i) + (dts/dxs^2)*Ds*(Z(j,i+1) - 2*Z(j,i) + Z(j,i-1))-((Vmax*dts*Z(j,i))/(Km+Z(j,i))); end Z(j+1,1)=Z(j,1)+dts*Ds*2*(Z(j,2)-Z(j,1))./dxs.^2-((Vmax*dts*Z(j,i))/(Km+Z(j,i))); %Neumann Boundary Condition %2nd Derivative Central Difference Product Iteration% for i=2:numx-1 Z(j+1,i) = Z(j,i) + (dtp/dxp^2)*Dp*(Z(j,i+1) - 2*Z(j,i) + Z(j,i-1))+((Vmax*dtp*Z(j,i))/(K