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

I am trying to integrate the following set of PDE

Equation 8 presents a 4th order spacial derivative, so my u vector is u = [C,rho,phi,P] where P is the second derivative of phi with respect to x. P = phi/xx is the fourth equation of the PDE system. Equations 9,10, 11 and 12 are boundary conditions valid for both boundaries, x = +1, x = -1. The initial condition for P is assumed to be P(0,t) = 0.

function [c,f,s] = pdefun_bsk(xmesh,tspan,u,dudx)
eps = 0.01;
nu = 0.5;
d_c = 10;
l_c = eps*sqrt(d_c);
c = [1; 1; 0; 0];
f = [eps/d_c*(dudx(1) + u(2)*dudx(3) + nu*u(1)/(1-nu*u(1))*dudx(1)); eps/d_c*(dudx(2) + u(1)*dudx(3) + nu*u(2)/(1-nu*u(1))*dudx(1)); dudx(3); dudx(4)];
s = [0; 0; -u(4);-u(2)/eps^2-u(4)/l_c^2];
end

function [pl,ql,pr,qr] = bcfun_bsk(xl,ul,xr,ur,tspan)
v = 0.1;
pl = [0; 0; ul(3)+v; 0];
ql = [1; 1; 0; 1]; 
pr = [0; 0; ur(3)-v; 0];
qr = [1; 1; 0; 1];
end

function u0 = icfun_bsk(xmesh)
v = 0.1;
u0 = [1; 0; v*xmesh; 0];
end

clear
clc
%% Data
m = 0;
xmesh = linspace(-1,1,1000);
tend = 2;
tspan = linspace(0,tend,1000);
tlength = length(tspan);
xlength = length(xmesh);
%% Solution
SOL = pdepe(m,@pdefun_bsk,@icfun_bsk,@bcfun_bsk,xmesh,tspan);

I get the following error

Warning: Failure at t=5.721307e-01.  Unable to meet integration tolerances without reducing the step
size below the smallest value allowed (1.776357e-15) at time t. 
> In ode15s (line 668)
  In pdepe (line 289) 
Warning: Time integration has failed. Solution is available at requested time points up to
t=5.705706e-01. 
> In pdepe (line 303) 

Is the problem due to the fact that I am trying to solve a 4-th order PDE? I tried to increase the tolerance (AbsTol e RelTol) to a few units and I do get a solution without errors, but it seems to be wrong.

[Bill Greene]

If you integrate to tend=.56 and then plot your fourth dependent variable at the right end as a function of time, you will see that it is rapidly approaching minus infinity.

[Alessandro Ferreri]

Hi Bill,

Thank you very much for the quick reply! You are right, the fourth variable does diverge on the right-hand side. I found this set of equations in a paper and a solution is presented. The set was solved with a different software though. Do you think the problem is numerical, or is there an issue with the equations themselves? I also tried to increase the number of steps in both time and space vectors up to 10000, without any luck.

Answer 1

[Bill Greene]

I have occasionally seen similar problems with pdepe in the past.

I have written a PDE solver that has input similar to pdepe but

does not exhibit this spurious divergence. If you want to try it,

you can download it here, pde1dM, and simply change your call from pdepe to pde1dM.

As an aside, the number of points in tspan has no effect on the solution

accuracy so I suggest you use far less than the 1000 you have chosen.

The number of points in xmesh does affect the accuracy but 1000 is

a very fine mesh; I suggest trying say, 50, and then increase that to

say, 100 and see how it affects the solution.

[Alessandro Ferreri]

Hi Bill,

thank you for the suggestion about xmesh and tspan vectors.

I tried your pde1dM function and in fact it yields a solution without any errors. It is the same solution that I get if I increase the relative tolerance in pdepe up to 30. However, this behavior differs from what I find in the original paper. I checked the way the equations are written and the parameters many times.

The right solution should be the one presented in the following paper in fig 7. EDIT: The most noticeable difference is that the third variable, phi, changes significantly between t = 0 and t = 2, whereas in my solution it barely changes.

Zhao, H., 2011. Diffuse-charge dynamics of ionic liquids in electrochemical systems. Physical Review E, 84(5), p.051504.