PDEPE: Unable to meet integration tolerances without reducing the step size below the smallest value allowed at time t.

Hello,

I am working on solving a system of 3 PDEs using the pdepe solver. I am receiving an error message that it is "unable to meet integration tolerances" which from what I have read, could mean a number of things may be wrong with my setup for using pdepe. I have tried altering a few different things but with no success.

The .pdf above details the system of PDEs that I wish to solve along with the initial and boundary conditions.

One place that I believe may be causing the error is in the boundary conditions in my j term. The term calls for the n variable evaluate at x = 0 but I was unsure if using ul(1) was the correct way to implement that. When I removed that term altogether, I still received a warning saying "unable to meet integration tolerances" but that it failed at a value of t that was slightly greater than the value of t it failed at whilst the ul(1) term was still in place.

Any help or insight is greatly appreciated!

MATLAB code:

global d alpha phi D_0 D_I D_I3 k_edr k_eer k_rg k_ext I2_0 I_0 Keq b gamma theta beta m1 n_0 I3_0 C1 C2 C3 e0 

e0 = 1.602e-19; % elementary charge (Coulombs)
d = 12.5e-4;    % cm
alpha = 0.2e-4; % cm^-1
phi = 6.1e16;   % cm^-2 s^-1
D_0 = 0.4;      % cm^2 s^-1
D_I = 4.4e-6;   % cm^2 s^-1
D_I3 = 3.6e-6;  % cm^2 s^-1
k_edr = 3e-11;  % cm^3 s^-1
k_eer = 2e-16;  % cm^3 s^-1
k_rg = 1.3e-15; % cm^3 s^-1
k_ext = 10^12;  % cm s^-1
n_0 = 1e14;     % cm^-3
I2_0 = 0.05;    % M
I_0 = 0.6;      % M
Keq = 10^8.3;
b = 1;           
gamma = 1.4;    % gamma/omega^2  
theta = 0.45;
beta = 0.35;
m1 = 1.55;

% pde coefficients for ease of use
C1 = gamma*(1-theta)*D_0;
C2 = gamma*(theta/1.5)*D_I;
C3 = gamma*(theta/1.5)*D_I3;

% solving for the equilibrium concentrations of I3_0 and I_0
syms I3
eqn = solve(I3/((I_0 - I3)*(I2_0 - I3)) == Keq, I3);
I3_0 = double(eqn(1));
I_0 = I_0 - I3_0;


x = linspace(0,d,101);
t = linspace(0,10,10);
m = 0;

sol = pdepe(m,@pde,@pde_ic,@pde_bc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3 = sol(:,:,3);


%---Functions----------------------------------------
% u(1) = n, u(2) = [I-], u(3) = [I3-]

function [c,f,s] = pde(x,t,u,dudx)
global phi alpha k_eer k_edr k_rg C1 C2 C3
S = phi.*alpha.*exp(-alpha.*x)./(k_rg.*u(2)+k_edr.*u(1));

c = [0; 0; 0];
f = [C1; C2; C3] .* dudx;
s = [phi.*alpha.*exp(-alpha.*x) - k_eer.*u(3).*u(1) - k_edr.*S.*u(1)
    1.5.*k_eer.*u(3).*u(1) - 1.5.*k_rg.*S.*u(2)
    -0.5.*k_eer.*u(3).*u(1) - 0.5.*k_rg.*S.*u(2)];
end

function u0 = pde_ic(x)
global n_0 I_0 I3_0
u0 = [n_0; I_0; I3_0];
end

function [pl,ql,pr,qr] = pde_bc(xl,ul,xr,ur,t)
global C1 C2 C3 e0 theta D_0 k_ext gamma
j = e0.*k_ext.*ul(1);

pl = [-j./(e0.*(1-theta).*D_0.*gamma); 0; 0];
ql = [1./C1; 1; 1];
pr = [0; 1.5.*j./(e0.*(1-theta).*D_0.*gamma); -0.5.*j./(e0.*(1-theta).*D_0.*gamma)];
qr = [1; 1./C2; 1./C3];
end

Warning messages:

When ul(1) is used:

Warning: Failure at t=2.357029e-02.  Unable to meet integration tolerances without reducing the step size below the smallest value
allowed (5.551115e-17) at time t. 
> In ode15s (line 653)
  In pdepe (line 289)
  In Fig3 (line 43) 
Warning: Time integration has failed. Solution is available at requested time points up to t=0.000000e+00. 
> In pdepe (line 303)
  In Fig3 (line 43) 

When ul(1) is removed:

Warning: Failure at t=3.909675e-02.  Unable to meet integration tolerances without reducing the step size below the smallest
value allowed (1.110223e-16) at time t. 
> In ode15s (line 653)
  In pdepe (line 289)
  In Fig3 (line 43) 
Warning: Time integration has failed. Solution is available at requested time points up to t=0.000000e+00. 
> In pdepe (line 303)
  In Fig3 (line 43)  

Answer 1

[Bill Greene]

You are trying to solve a system of ODE which is something pdepe is not designed for. Specifically, the pdepe documentation says, "At least one equation must be parabolic.". This means that at least one of your c-coefficients must be non-zero. I'm surprised pdepe does not specifically test for this and issue an error message.

I suggest you look at one of the boundary value solver functions, instead (e.g. bvp4c).

[Wes Fermanich]

Thank you very much for this answer. It wasn't apparent to me from the error message that I was getting that pdepe wasn't designed to handle this kind of steady state system. This clears it up!

[Bill Greene]

The documentation is not very clear on this point and the error message doesn't help at all.