Error using pdepe - Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative.

Hello,

I'm using pdepe to solve for a reaction-diffusion equation. When I set the initial concentration to anything between 0-250, it runs fine, but the moments, I get this error saying

"Error using pdepe Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative."

Am I setting my boundary conditions incorrectly? I want the left hand side (at x = 0, C_L = C_L0) and on the right hand side I want the flux to equal to 0 (at x = x_f, dC_L/dx = 0).

function time_dependent_pdepe
clear all; close all; clc;

D_ij = 30;
L0 = 1000; %c0 [nM]
x_f =40; %Length of domain [um]
maxt = 2; %Max simulation time [s]

m = 0; 
x = linspace(0,x_f,1000); %xmesh
t = linspace(0,maxt,50) %tspan

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

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

end
title('Time Dependent (Diffusion Only)')
xlabel('Distance (\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*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 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));

% PDE
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)
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

Where am I going wrong?

Thank you.

[Bill Greene]

Your post doesn't say how you reproduce this error?

[Jennifer Yang]

I get this error whenever I set L0 (line 5) to anything higher than 250.

Answer 1

[Bill Greene]

The fundamental cause of your problem is that your initial conditions are inconsistent with your boundary condition at x=0. Mathematically this is referred to as an ill-posed problem and the PDE has no solution.

Sometimes, if the inconsistency is not too large, the numerical algorithm in pdepe can adjust the initial conditions to make them consistent with the boundary conditions. But, in your case, for large values of L0, it is not able to do this.

One simple way to fix this problem is to modify your initial condition function slightly:

function u0 = DiffusionICfun(x)
  if(x==0)
    u0 = L0;
  else
     u0 = 0;
  end
end