PDEPE: Unable to meet integration tolerances without reducing the step size below the smallest value allowed at time t.
3 views (last 30 days)
Show older comments
Wes Fermanich
on 24 Jun 2020
Commented: Bill Greene
on 25 Jun 2020
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)
0 Comments
Accepted Answer
Bill Greene
on 25 Jun 2020
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).
2 Comments
Bill Greene
on 25 Jun 2020
The documentation is not very clear on this point and the error message doesn't help at all.
More Answers (0)
See Also
Categories
Find more on Eigenvalue Problems in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!