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

13 visualizzazioni (ultimi 30 giorni)
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.
  2 Commenti
Bill Greene
Bill Greene il 13 Mar 2021
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
Alessandro Ferreri il 14 Mar 2021
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.

Accedi per commentare.

Risposte (1)

Bill Greene
Bill Greene il 14 Mar 2021
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.
  1 Commento
Alessandro Ferreri
Alessandro Ferreri il 15 Mar 2021
Modificato: Alessandro Ferreri il 15 Mar 2021
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.

Accedi per commentare.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by