# Solving PDE with unreasonable results

I have a simplified version of my chemical reaction problem. Basically it solves the concentration (ca) and grain (g2) diameter. The ca function is depended on particle diameter R, while g2 function is depend on time t. The ca and g2 function are depended on each others. In this case, both ca and g2 should be a function of R and t. For my understanding, it should be clarified as PDE problem and can be solved using MethodofLines.

The code is given as:

del = 10^(-8);

Monitor[{casol, gsol} =

NDSolveValue[   {
D[R^2*D[ca[R, t], R], R]/R^2 -
( 39.7*(g2[R, t])^2*ca[R, t])/(1 + g2[R, t]*(1-g2[R, t])) == 0,

D[g2[R, t], t] == -ca[R, t]/(1 + g2[R, t]*(1 - g2[R, t])),

ca[R, 0.0] == 1.0,
ca[1.0, t] == 1.0,
(D[ca[R, t], R] /. R -> del) == 0,

g2[R, 0.0] == 1.0  },

{ca, g2}, {R, del, 1}, {t, 0, 1} ,

Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
MinPoints -> 200}} ,

EvaluationMonitor :> (monitor =
Row[{"t=", CForm[t], " gsol=", CForm[g2[1.0, t]]}])],

monitor]


It can solve the problem with a few warnings but the result is not right. The g2 gives you a constant value of 1 with (R,t). This is not right. A simple analytic analysis for g2 equation at R=1 can prove g2 is a function of t, and should not be constant.

Also it seems the boundary/ initial conditions are not forced in the results. For example, if I try to Plot[casol[1.0, t], {t, 0, 1}]. It expect give you value of 1,as this is part of the boundary condition, but the result shows some unreasonable value.

Not sure what cause the problems? Really appreciate some one can provide some hints.

• "Warning: an insufficient number of boundary conditions have been specified for the direction of independent variable R. Artificial boundary effects may be present in the solution." That's a warning that you should not take light hearted. – Henrik Schumacher Jul 31 '19 at 4:02
• Moreover, I don't see any timederivatives of ca. Are you sure that your entered the PDE correctly? – Henrik Schumacher Jul 31 '19 at 4:04
• Yes,I initial suspected this as well. But from the chemical reaction point of view, this is called pseudo-steady state. The mathematics is taken from Chemical Engineering Science, Vol 34, 1979,1072-1075 with some simplification. This make me thought this type of problem should not use MethodOfLines to solve, but have no idea what other method can be used. – Kim Tang Jul 31 '19 at 4:14

The system of equations can be solved using the MethodOfLines, but several options need to be added.The physical model is not finalized, since in the process of calculation for t>= 0.44 we have g2<0. You can calculate up to t = 0.49 after which an instability arises.

del = 10^(-8); tmax = 0.49;

Monitor[{casol, gsol} =
NDSolveValue[{D[R^2*D[ca[R, t], R], R]/
R^2 - (397/10*(g2[R, t])^2*ca[R, t])/(1 +
g2[R, t]*(1 - g2[R, t])) == 0,
D[g2[R, t], t] == -ca[R, t]/(1 + g2[R, t]*(1 - g2[R, t])),
ca[R, 0] == 1, ca[1, t] == 1, Derivative[1, 0][ca][del, t] == 0,
g2[R, 0] == 1}, {ca, g2}, {R, del, 1}, {t, 0, tmax},
Method -> {"IndexReduction" -> Automatic,
"EquationSimplification" -> "Residual",
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid",
"MinPoints" -> 137, "MaxPoints" -> 137,
"DifferenceOrder" -> "Pseudospectral"}}},
EvaluationMonitor :> (monitor =
Row[{"t=", CForm[t], " gsol=", CForm[g2[1.0, t]]}])], monitor]

{Plot3D[casol[x, t], {x, del, 1}, {t, 0, .49}, PlotRange -> All,
AxesLabel -> Automatic, Mesh -> None, ColorFunction -> "Rainbow"],
Plot3D[gsol[x, t], {x, del, 1}, {t, 0, .49}, PlotRange -> All,
AxesLabel -> Automatic, Mesh -> None, ColorFunction -> "Rainbow"]}


• Alex: Thanks very much. I know my another question is hard to answer, but it is very useful for a learner and will be much appreciated. How to get to all this options to get the case running? are there any tutorial or reference that you can recommend? It seems we need to dig the details of the algorithm, which I do not mind but it is obvious against the principle of Wolfram language. Thanks again. – Kim Tang Jul 31 '19 at 7:03
• @KimTang I found these options to solve another problem. But in this case also works. There is an example in the documentation, but I cannot find it now. – Alex Trounev Jul 31 '19 at 7:38