# NDSolveValue and pdepe of MATLAB disagree under spherical coordinates

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=du*8.3E-10;
s=-0.01*u;
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}]

• For reference, could you include the MATLAB code you were using (if it isn't too long)? Commented Feb 11, 2021 at 6:09
• @J.M. I have post my Matlab code. Thanks! Commented Feb 11, 2021 at 7:44
• Thanks for your solutions. Commented Feb 15, 2021 at 2:23
• As mentioned in my last comment, you should modify the inner b.c. to Derivative[0,1][CC][t,ri/50]==0 to avoid the removable singularity. With this b.c., the output is consistent with that of MATLAB. Also, why do you type {r,0,ri,ri/50}? This syntax isn't documented and colored red in the notebook. Though NDSolve still (surprisingly) gives a result, the mysterious syntax clearly causes certain problem. Commented Feb 21, 2021 at 9:11
• It should be {r, ri/50, ri}. Commented Feb 21, 2021 at 12:10

## 1 Answer

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:

• Good find! Did not know this existed. Commented Feb 22, 2021 at 10:39
• I failed to find the "Coordinates" option in Documentation Center, which might be the reason that I did not follow your point. Do you mean that "Coordinates" option is used to make a custom spatial grid ? Commented Feb 23, 2021 at 11:51
• @hithere It's discussed in the tutorial The Numerical Method of Lines. Yes, it's used for customizing spatial grid. Commented Feb 23, 2021 at 12:03