NDSolveValue and pdepe of MATLAB disagree under spherical coordinates

Question

Assuming there is a spherical activated carbon particle, the concentration outside the grain is 200 (C[t,ri] == 200), that inside the grain is 0 (Derivative[0, 1][C][t, 0] == 0), the initial condition for PDE is C[0, r] == 0, the radius for the particle is ri(0.001 m). I am new to Mathematica, and when I run this code, the plot shows slight difference with that in MATLAB.

The plot in Mathematica sees a decline of concentration resulting a negative concentration after a period of adsorption, while it remains steady in MATLAB figure. Actually after adsorption equilibrium, the concentration should be stable . Is there a way to fix my code? Thanks.

The concentration becomes negative with growing time:

Figure in MATLAB:

The equation is like this:

Mathematica code:

TT = 3000; ri = 0.001;
sol = NDSolveValue[{D[CC[t, r], t] == 
    8*10^-10/r^2*D[r^2*D[CC[t, r], r], r] - 0.01*CC[t, r], CC[0, r] == 0, 
  CC[t, ri] == 200, Derivative[0, 1][CC][t, 0] == 0}, CC, {t, 0, TT}, {r, 0, ri, ri/50}]
Table[Plot3D[sol[t, x], {x, 0, ri}, {t, 0, VT}, PlotRange -> Full, 
  ViewPoint -> {-5, -2, 3}], {VT, 100, TT, TT/5}]

This is my MATLAB code:

function test
global c0;
c0=200;         
tc=1200;
ri=0.001;
dx=70;
x=linspace(0,ri,dx); 
dtL=25;dtR=10;

t=[linspace(0,tc/4,dtL) linspace(tc/4+1,2000000,dtR)];
m=2;
sol=pdepe(m,@mpde,@mpic,@mpbc,x,t);
c=sol(:,:,1);
figure(1);
fig1=surf(x,t,c);
xlabel(' x');
ylabel(' t');
end
function [c,f,s]=mpde(x,t,u,du) 
c=1;
f=du8.3E-10;
s=-0.01u;
end
function [pl,ql,pr,qr]=mpbc(xl,ul,xr,ur,t) 
global c0;
pl=0;
ql=1;
pr=ur-c0;
qr=0;
end
function u0=mpic(x)
u0=0;
end

Update 2.21:

Some results actually confused me that,

1. the stable value for r = 0 in MATLAB is around 43, while it is 22 in Mathematica.

2. a gradually increasing could be seen in initial 100 seconds, while it instantly reaches stable in Mathematica.

Does this discrepancy result from the difference between MATLAB and Mathematica or the code itself?

T2 = 1000;
sol = NDSolveValue[{D[CC[r, t], t] == 
    8.3 10^-10/r^2*D[r^2*D[CC[r, t], r], r] - 0.01 CC[r, t], CC[r, 0] == 0, 
   CC[ri, t] == 200, Derivative[1, 0][CC][0, t] == 0},CC, {t, 0, T2}, {r, 0, ri, ri/50},
   Method -> {"MethodOfLines", 
    "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 100}}]
Table[Plot3D[sol[x, t], {x, 0, ri}, {t, 0, VT2}, PlotRange -> All, 
  ColorFunction -> ColorData["BlueGreenYellow"], PlotPoints -> 100, 
  PlotLegends -> BarLegend[Automatic, LegendMarkerSize -> {10, 200}], TicksStyle -> 13, 
  ViewPoint -> {-3, -2, 1}, AxesLabel -> {"r", "t", "z"}], {VT2, 100, T2, T2/3}]
ContourPlot[sol[r, t], {r, 0, ri}, {t, 0, T2}]

Answer 1

Interesting. OP has hit on an undocumented syntax of NDSolve that seems never be discussed in this site.

As we can see, OP has typed {r, 0, ri, ri/50} inside NDSolve. This syntax isn't mentioned in the document AFAIK, and is colored red at least in the notebook of v8 and v9, but surprisingly, NDSolve still gives an output:

T2 = 1000; ri = 0.001;
inner = 0;
solaccident = 
 NDSolve[{D[CC[r, t], t] == 8.3 10^-10/r^2*D[r^2*D[CC[r, t], r], r] - 0.01 CC[r, t], 
     CC[r, 0] == 0, CC[ri, t] == 200, Derivative[1, 0][CC][inner, t] == 0}, 
    CC, {t, 0, T2}, {r, inner, ri, ri/50}][[1, 1, -1]] // Quiet

How does NDSolve interprete the {r, 0, ri, ri/50}? It turns out to be equivalent to "Coordinates" -> {inner, ri, ri/50}!:

solcheck = NDSolve[{D[CC[r, t], t] == 
     8.3 10^-10/r^2*D[r^2*D[CC[r, t], r], r] - 0.01 CC[r, t], CC[r, 0] == 0, 
    CC[ri, t] == 200, Derivative[1, 0][CC][inner, t] == 0}, 
   CC, {t, 0, T2}, {r, inner, ri}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "Coordinates" -> {inner, ri, ri/50}}}][[1, 1, -1]]
solcheck == solaccident
(* True *)

Clearly, the spatial grid {0, ri/50, ri} is way too coarse, and results in an inaccurate result.

Remark

• The option "Coordinates" is documented in the tutorial The Numerical Method of Lines, particularly Controlling the Spatial Grid Selection subsection.

The following is the simplest way to fix the code:

T2 = 1000; ri = 0.001;
inner = ri/50;
soltraditional = 
 NDSolveValue[{D[CC[r, t], t] ==8.3 10^-10/r^2*D[r^2*D[CC[r, t], r], r] - 0.01 CC[r, t],
    CC[r, 0] == 0, CC[ri, t] == 200, Derivative[1, 0][CC][inner, t] == 0}, 
  CC, {t, 0, T2}, {r, inner, ri}]

Though ibcinc warning pops up, it's not causing problem here, because the Neumann b.c. isn't the inconsistent one. (For more info check this post. )

I've moved the inner boundary to ri/50 to avoid the removable singularity at origin. In principle this isn't necessary, but NDSolve will return the input otherwise, which may be a bug. Is it possible to circumvent the bug without moving the inner boundary? Once again, it's surprisingly possible! We just need to set the "Coordinates" option. Let's make use of the hidden syntax we just found:

T2 = 1000; ri = 0.001;
inner = 0; dr = ri/25;
soltest = NDSolveValue[{D[CC[r, t], t] == 
    8.3 10^-10/r^2*D[r^2*D[CC[r, t], r], r] - 0.01 CC[r, t], CC[r, 0] == 0, 
   CC[ri, t] == 200, Derivative[1, 0][CC][inner, t] == 0}, CC, {t, 0, T2}, 
  Flatten@{r, Range[inner, ri, dr]}]

Power::infy and Infinity::indet will pop up, but don't worry.

---

Further check shows this hidden syntax is introduced in v3: