5
$\begingroup$

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:

enter image description here

Figure in MATLAB:

enter image description here

The equation is like this:

enter image description here

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}]

enter image description here enter image description here

$\endgroup$
8
  • $\begingroup$ For reference, could you include the MATLAB code you were using (if it isn't too long)? $\endgroup$ Commented Feb 11, 2021 at 6:09
  • $\begingroup$ @J.M. I have post my Matlab code. Thanks! $\endgroup$
    – hithere
    Commented Feb 11, 2021 at 7:44
  • $\begingroup$ Thanks for your solutions. $\endgroup$
    – hithere
    Commented Feb 15, 2021 at 2:23
  • 1
    $\begingroup$ 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. $\endgroup$
    – xzczd
    Commented Feb 21, 2021 at 9:11
  • 1
    $\begingroup$ It should be {r, ri/50, ri}. $\endgroup$
    – xzczd
    Commented Feb 21, 2021 at 12:10

1 Answer 1

6
$\begingroup$

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

enter image description here

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:

enter image description here

$\endgroup$
3
  • $\begingroup$ Good find! Did not know this existed. $\endgroup$
    – user21
    Commented Feb 22, 2021 at 10:39
  • $\begingroup$ 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 ? $\endgroup$
    – hithere
    Commented Feb 23, 2021 at 11:51
  • 2
    $\begingroup$ @hithere It's discussed in the tutorial The Numerical Method of Lines. Yes, it's used for customizing spatial grid. $\endgroup$
    – xzczd
    Commented Feb 23, 2021 at 12:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.