trouble using pdepe to solve system of pdes

Hello,

I have to solve the following system of pdes:

The code below uses pdepe to solve it, but it returned the error:

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

Could anyone help solve the problem?

Thank you in advance!

global D1 H D2
H = 91e-6; % paper height
D1 = 3e-6; % diffusion constant of pores
D2 = 1e-14;
t0 = 20

t0 = 20

global T k M0 C K eta rho cw cf0
T = 298; 
k = 0.0035;
M0 = 0.0329; 
C = 39.09; 
K = 0.865; 
eta = 0.47; 
rho = 1500; 
cf0 = 0.25*rho; 
cw = 1000; 

global csp
A1 = 1.3258;
A2 = -0.003931;
A3 = 20.7115;
A4 = -5364.05;
A5 = -17.58;
csp = vpa((A1/T)*exp((A2*T^2 + A3*T + A4)/(T + A5)));

tSpan1 = linspace(0,t0,1001);
xmesh = linspace(0,H,20);
m = 0;
sol1 = pdepe(m,@pdefun,@icfun,@bcfun1,xmesh,tSpan1);

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

function [c,f,s] = pdefun(x,t,u,dudx)
global eta D1 rho M0 C K csp k D2
c = [1-eta;eta];
f = [(1-eta)*D2;eta*D1].*dudx;
rh = u(2)/csp;
A = rho*M0*C*K*(1/csp)/((1-K*rh)*(1-K*rh+C*K*rh));
s = k*[A*u(2) - u(1);-A*u(2) + u(1)];
end

function u0 = icfun(x)
u0 = [0;0];
end

function [pl,ql,pr,qr] = bcfun1(xl,ul,xr,ur,t)
global cw
pl = [ul(1);ul(2)-cw];
ql = [0;0];
pr = [0;0];
qr = [1;1];
end

Answer 1

[Torsten]

Your initial condition for u(2) at x = 0 is not consistent with your boundary condition.

Use

function u0 = icfun(x)
  global H cw
  u0 = [0;0];
  if x==0
    u0(2)=cw;
  end
end

instead.

And remove the vpa in the evaluation of csp.

[Chin Ching]

Hi Torsten,

I just made the adjustments. However the same error still popped up.

[Torsten]

For me it works in the current online MATLAB release R2024b. What MATLAB version do you use ?

global D1 H D2
H = 91e-6; % paper height
D1 = 3e-6; % diffusion constant of pores
D2 = 1e-14;
t0 = 20;

global T k M0 C K eta rho cw cf0
T = 298; 
k = 0.0035;
M0 = 0.0329; 
C = 39.09; 
K = 0.865; 
eta = 0.47; 
rho = 1500; 
cf0 = 0.25*rho; 
cw = 1000; 
global csp
A1 = 1.3258;
A2 = -0.003931;
A3 = 20.7115;
A4 = -5364.05;
A5 = -17.58;
csp = (A1/T)*exp((A2*T^2 + A3*T + A4)/(T + A5));
tSpan1 = linspace(0,t0,1001);
xmesh = linspace(0,H,20);
m = 0;
sol1 = pdepe(m,@pdefun,@icfun,@bcfun1,xmesh,tSpan1);
u1 = sol1(:,:,1);
u2 = sol1(:,:,2);
plot(tSpan1,u2(:,end))

function [c,f,s] = pdefun(x,t,u,dudx)
global eta D1 rho M0 C K csp k D2
c = [1-eta;eta];
f = [(1-eta)*D2;eta*D1].*dudx;
rh = u(2)/csp;
A = rho*M0*C*K*(1/csp)/((1-K*rh)*(1-K*rh+C*K*rh));
s = k*[A*u(2) - u(1);-A*u(2) + u(1)];
end
function u0 = icfun(x)
  global cw
  u0 = [0;0];
  if x==0
    u0(2)=cw;
  end
end
function [pl,ql,pr,qr] = bcfun1(xl,ul,xr,ur,t)
global cw
pl = [ul(1);ul(2)-cw];
ql = [0;0];
pr = [0;0];
qr = [1;1];
end

[Chin Ching]

it works for me now. Thanks a lot Torsten!