0
$\begingroup$

I'm trying to use MATLAB solve a reaction-diffusion equation in both steady state (using bvp4c) and time dependency (using pdepe). When plotting both of the solutions, I do not get the same answer. For my boundary conditions I have: at x = 0, C = 0 and at x = L, dC/dx = 0.

For my steady state case (bvp4c):

  function steady_state_bvp4c
  close all; clear all; clc;

  %Diffusion coefficient
  D_ij = 30; %3.0*10^-7 cm^2/s -> 30 um^2/s

  %Initial concentration at x=0
  L0 = 50;

  %Total length
  x_f = 20;

  %Rate constants
  k_1 = 0.00193;
  k_2 = 0.00255;
  k_3 = 4.09;
  d_1 = 0.007;
  d_2 = 3.95*10^-5;
  d_3 = 2.26;

  %Equilibrium constants
  K_1 = 3.63;
  K_2 = 0.0155;
  K_3 = 0.553;
  K_4 = 9.01;
  K_i = 0.139;

  %Total Receptors
  N_T =1.7;

  solinit = bvpinit(linspace(0,x_f,11),[0.5 0]);

  sol = bvp4c(@(x,y)odefcn(x,y,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T),@twobc,solinit);

  figure(1)
  plot(sol.x,(sol.y(1,:)),'LineWidth',1)
  title('Steady State')
  xlabel('Distance from Device (\mum)')
  ylabel('Concentration (nM)')
  axis([0 x_f 0 L0])

  function dy = odefcn(x,y,D_ij,k_1,k_2,k_3,d_1,d_2,d_3,K_1,K_2,K_3,K_4,K_i,N_T)
  C_L    = y(1);
  dC_Ldx = y(2);

  N_r = -(K_3.*K_4.*K_i.*(K_2 + C_L - ((8.*K_2.*C_L.^2.*N_T + 8.*K_2.*K_3.*C_L.*N_T + K_2.^2.*K_3.*K_4.*K_i + K_3.*K_4.*K_i.*C_L.^2 + 8.*K_2.^2.*K_3.*K_i.*N_T + 2.*K_2.*K_3.*K_4.*K_i.*C_L + 8.*K_2.*K_3.*K_i.*C_L.*N_T)./(K_3.*K_4.*K_i)).^(1./2)))./(4.*C_L.^2 + 4.*K_3.*C_L + 4.*K_2.*K_3.*K_i + 4.*K_3.*K_i.*C_L);
  N_R = (C_L./K_1).*N_r;
  N_r2 = (1./K_4).*(N_r).^2;
  N_rR = (C_L./(2.*K_2.*K_4)).*(N_r).^2;
  N_pR = (C_L./(2.*K_2.*K_4.*K_i)).*(N_r).^2;
  N_R2 = ((C_L.^2)./(K_2.*K_3.*K_4.*K_i)).*(N_r).^2;

  R_L = ((2.*d_1.*(N_T-N_r-N_r2-N_rR-N_pR-N_R2))-(2.*k_1.*C_L.*N_r))+...
      ((2.*d_2.*(N_T-N_r-N_R-N_r2-N_pR-N_R2)))-((k_2.*C_L.*(N_T-N_r-N_R-N_rR-N_pR-N_R2)))+...
      (d_3.*(N_T-N_r-N_R-N_r2-N_rR-N_pR))-(2.*k_3.*C_L.*(N_T-N_r-N_R-N_r2-N_rR-N_R2));

  dy = zeros(2,1);
  dy(1) = dC_Ldx;
  dy(2) = (-R_L/D_ij);
  end

  function res = twobc(ya,yb)
  res = [ ya(1)-(L0); yb(2) ];
  end
  end

steady state plot

For my time dependent case:

  function time_dependent_pdepe
  clear all; close all; clc;

  D_ij = 30; %Diffusion coefficient D (3.0*10^-7 cm^2/s -> 30 um^2/s)
  L0 = 50; %c0 [nM]
  x_f =20; %Length of domain [um]
  maxt = 1; %Max simulation time [s]

  m = 0; %Parameter corresponding to the symmetry of the problem
  x = linspace(0,x_f,100); %xmesh
  t = linspace(0,maxt,100); %tspan

  sol = pdepe(m,@DiffusionPDEfun,@DiffusionICfun,@DiffusionBCfun,x,t,[]);
  u = sol;

  % Plotting
  hold all
  for n = linspace(1,length(t),10)
  plot(x,sol(n,:),'LineWidth',2)

  end
  title('Time Dependent')
  xlabel('Distance from Device (\mum)')
  ylabel('Concentration (nM)')
  axis([0 x_f 0 L0])



  function [c,f,s] = DiffusionPDEfun(x,t,u,dudx)

  D = D_ij;

  %Rate constants
  k_1 = 0.00193;
  k_2 = 0.00255;
  k_3 = 4.09;
  d_1 = 0.007;
  d_2 = 3.95;
  d_3 = 2.26;

  %Equilibrium constants
  K_1 = 3.63;
  K_2 = 0.0155;
  K_3 = 0.553;
  K_4 = 9.01;
  K_i = 0.139;

  %Total Recptors
  N_T = 1.7;

  N_r = -(K_3.*K_4.*K_i.*(K_2 + u - ((8.*K_2.*u.^2.*N_T + 8.*K_2.*K_3.*u.*N_T + K_2.^2.*K_3.*K_4.*K_i + K_3.*K_4.*K_i.*u.^2 + 8.*K_2.^2.*K_3.*K_i.*N_T + 2.*K_2.*K_3.*K_4.*K_i.*u + 8.*K_2.*K_3.*K_i.*u.*N_T)./(K_3.*K_4.*K_i)).^(1./2)))./(4.*u.^2 + 4.*K_3.*u + 4.*K_2.*K_3.*K_i + 4.*K_3.*K_i.*u);
  N_R = (u./K_1).*N_r;
  N_r2 = (1./K_4).*(N_r).^2;
  N_rR = (u./(2.*K_2.*K_4)).*(N_r).^2;
  N_pR = (u./(2.*K_2.*K_4.*K_i)).*(N_r).^2;
  N_R2 = ((u.^2)./(K_2.*K_3.*K_4.*K_i)).*(N_r).^2;

  R_L = ((2.*d_1.*(N_T-N_r-N_r2-N_rR-N_pR-N_R2))-(2.*k_1.*u.*N_r))+...
      ((2.*d_2.*(N_T-N_r-N_R-N_r2-N_pR-N_R2)))-((k_2.*u.*(N_T-N_r-N_R-N_rR-N_pR-N_R2)))+...
      (d_3.*(N_T-N_r-N_R-N_r2-N_rR-N_pR))-(2.*k_3.*u.*(N_T-N_r-N_R-N_r2-N_rR-N_R2));

  c = 1;
  f = D_ij.*dudx;
  s = R_L;
   end

  function u0 = DiffusionICfun(x)
  u0 = 0;
  end

  function [pl,ql,pr,qr] = DiffusionBCfun(xl,ul,xr,ur,t)
  % Boundary conditions for x = 0 and x = L;

  c0 = L0;

  % BCs: No flux boundary at the right boundary and constant concentration on
  % the left boundary

  % At x = 0, c0 = L0
  pl = ul-c0; 
  ql = 0; 

  % At x = L, dc/dx = 0
  pr = 0;
  qr = 1;

  end

  end

time dependent plot

Is there a difference in the way I am setting up both situations? I would greatly appreciate any help! Thank you!

$\endgroup$
  • $\begingroup$ There is way to much code for someone to go through, so you are unlikely to get an answer. I would suggest finding the absolute simplest setting where you get different answers for time-dependent and stationary code, and post that. Doing this type of exercise will probably illuminate the issue for you anyway. $\endgroup$ – Jeff Jul 17 '18 at 17:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.