# Problems with Numerical Solution - NDSolve

This is the code I use :

kd = 2;
n = 1.25;
T = 1000;
rw = 0.001;
R = rw + Sqrt[(T*Pi)/1.4];
a = 0.5;
b = 50;
c = 5.*10^-10;

sys = {(D[u[r, t], t] + a*(u[r, t] - m[r, t])) ==
(kd^(1/n)/n*D[u[r, t], {r, 2}] +
kd^(1/n)/r*D[u[r, t], r])*(-D[u[r, t], r])^((1 - n)/n),
b*D[m[r, t], t] == a*(u[r, t] - m[r, t]),
D[u[rw, t], r] == -(2^n/(rw^n*kd)),
u[R, t] == c,
u[r, c] == c,
m[r, c] == c};

sol = NDSolve[sys, {u, m}, {r, rw, R}, {t, c, T},
MaxStepSize -> {0.001, 1},
Method -> {"MethodOfLines",
"DifferentiateBoundaryConditions" -> {True,
"ScaleFactor" -> T}}]; // AbsoluteTiming

np =  LogLogPlot[u[1, t] /. sol, {t, 0.1, T}];
Show[ np,  PlotRange -> All]


But when the code is running, there are some messages:

NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent.
Power::infy: Infinite expression 1/0.^1.2 encountered.
Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.
General::stop: Further output of Power::infy will be suppressed during this calculation.


In the calculation, The memory is always insufficient and these messages only happened at n=1.25, 1.5, and 1.75. If n=1, the messages do not appear. However, The trend of all curves needs to be similar to n=1. So, how can I do to fix it?

• The problem is causing by the constant initial condtion, which is not allowed by terms like (-D[u[r, t], r])^((1 - n)/n). You need a non-constant initial condition, or modify the equation to avoid this singularity. – xzczd Nov 14 '17 at 11:12
• Can you provide some example, please?? I still cannot understand well about the constant initial condition. Do you mean that I need to change the eq., D[u[rw, t], r] == -(2^n/(rw^n*kd)) ?? – Hsin Nov 15 '17 at 2:25
• Notice that, since u[r, c] == c, D[u[r, c], r] == 0, so (-D[u[r, c], r])^((1 - n)/n) is no longer well defined when ((1 - n)/n) is a negative number. – xzczd Nov 15 '17 at 10:00