Pdepe: Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative.

Hi,

I was trying to solve PDEs system about diffusion on catalyst (modeled by deactivation model), which yields:

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

Error in Resolution (line 30)
sol = pdepe(m,@eqn2,@initial2,@bc2,x,t);

After read model in attchament, please see the code below:

global Def T Dp R cAO k kd kO cW alphaO gamma

%Parameters
T = 333; %[=]K
Dp = 2; %[=]m
Def = 10^-5; %
cAO = 0.0001; %[=]kmol/m^3 intial valure of concentration
cW = 0.0001;%[=]kmol/m^3
R = Dp/2;
k = (1.191*10^5)*exp(-7.544/(0.001607460438947*T));
kd = (4.299*10^3)*(exp(-6.86/(0.002102114315754*T)));
kO =k*cW;
gamma = 1;
alphaO = 1; %intial value of catalyst activity

%PDE2: MATLAB script M-file that solves the PDE 
%stored in eqn2.m, bc2.m, and initial2.m
m = 2; %spheric

x = linspace(0,R,10);
t = linspace(0,1,10);

sol = pdepe(m,@eqn2,@initial2,@bc2,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
%--------------------------------------------------------
%EQN2: MATLAB M-file that contains the coefficents for 
%a system of two PDE in time and one space dimension.
function [c,b,s] = eqn2(x,t,u,DuDx)

global Def kO kd gamma
c = [1; 1];
b = [Def; 0] .* DuDx;
s = [- kO*u(1)*(u(2));-kd*u(1)*((u(2)).^(gamma+1))];
% --------------------------------------------------------
%INITIAL2: MATLAB function M-file that defines initial conditions 
%for a system of two PDE in time and one space variable.
function value = initial2(x);

global cAO alphaO
value = [cAO; alphaO];
% --------------------------------------------------------
%BC2: MATLAB function M-file that defines boundary conditions 
%for a system of two PDE in time and one space dimension.
function [pl,ql,pr,qr] = bc2(xl,ul,xr,ur,t)

global cAO R
pl = [0; 0];
ql = [1; 0];
pr = [ur(R)-cAO; 0];
qr = [0; 0];

Answer 1

[Bill Greene]

The boundary conditions for your second PDE are invalid. These should work:

pl = [0; 0];
ql = [1; 1];
pr = [ur(1)-cAO; 0];
qr = [0; 1];

[Luigi Piero Di Bonito]

Thank you Sir, now the script work!

But I have another problem, I was trying to change Dp value from 2 [m] to 10^-6 [m], which yields:

Subscript indices must either be real positive integers or logicals.
Error in bc2 (line 10)
pr = [ur(Dp/2)-cAO; 0];

Error in pdepe (line 250)
[pL,qL,pR,qR] = feval(bc,xmesh(1),y0(:,1),xmesh(nx),y0(:,nx),t(1),varargin{:});

Error in Resolution (line 28)
sol = pdepe(m,@eqn2,@initial2,@bc2,x,t);

[Bill Greene]

Sorry, that is due to another error in your original script that I missed. I edited my answer to fix that.

[Luigi Piero Di Bonito]

Thank you very much Sir!