# solving a nonlinear PDE of infinite cylindrical domain with NDSolve

I'm trying to solve the following nonlinear PDE

Note that when n is equal to unit, the equation becomes linear. Thus, I applied the NDSolve function to solve it numerically. Here is my code:

  kd = 1; n = 2; T = 1000; R = 1 + Sqrt[T*Pi];
sys = {(-D[u[r, t], r])^(1/n - 1)*((kd^(1/n)/n)*D[u[r, t], r, r] + (kd^(1/n)/r)*D[u[r, t], r]) == D[u[r, t], t]
,
Derivative[1, 0][u][1, t] == -(2^n/kd), u[R, t] == 0, u[r, 0] == 0
};
AbsoluteTiming[sol = NDSolve[sys, {u}, {r, 1, R}, {t, 0, T},
Method -> {"MethodOfLines",
"DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> T}}] ]


The R value represent the influence radius to overcome the infinite boundary condition. It has been valid compared to reference solution (n=1). However, when I'm trying to evaluate the numerical result as n > 1 (in this case the n value is chosen as 2), the massages appear:

Power::infy: Infinite expression 1/0.^1.5 encountered. >>

Power::infy: Infinite expression 1/0.^1.5 encountered. >>

Power::infy: Infinite expression 1/0.^1.5 encountered. >>

General::stop: Further output of Power::infy will be suppressed during this calculation. >>

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

General::stop: Further output of Infinity::indet will be suppressed during this calculation. >>

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.. >>

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.. >>

NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.. >>

General::stop: Further output of NDSolve::ndnum will be suppressed during this calculation. >>


Which step I did wrong?

The reference solution is provided as the function

s[r, t]


The code is

  Vi[n_, i_] :=
Vi[n, i] = (-1)^(i + n/2) Sum[
k^(n/2) (2 k)! /( (n/2 - k)! k! (k - 1)! (i - k)! (2 k -
i)! ), { k, Floor[ (i + 1)/2 ], Min[ i, n/2] } ] // N;
Stehfest[F_, s_, t_, n_: 16] :=
If[n > 16, Message[Stehfest::optimalterms, n];
If[ OddQ[n], Message[Stehfest::odd, n];
"Enter an even number of terms",
If[n > 32, Message[Stehfest::terms, n];
" Try a smaller value for n. Maximum allowable n is 32 ",
Log[2]/t Sum[ Vi[n, i]*F /. s -> i Log[2]/t , {i, 1, n} ] ]],
If[ OddQ[n], Message[Stehfest::odd, n];
"Enter an even number of terms",
If[n > 32, Message[Stehfest::terms, n];
" Try a smaller value for n. Maximum allowable n is 32.",
Log[2]/t Sum[
Vi[n, i]*F /. s -> i Log[2]/t , {i, 1, n} ] ]]]  // N;
\[Beta]=0; rw=1;s[r_, t_] :=
Stehfest[(2^(2 + n) r^(1/2 - n/2)
rw BesselK[(-1 + n)/(-3 + n), -((
r^(3/2 - n/2) Sqrt[(
2^(1 + n) n p (\[Beta] + p \[Eta] + \[Beta] \[Eta]))/(
kd (\[Beta] + p \[Eta]))])/(-3 + n))])/(kd p (rw^2 Sqrt[(
2^(1 + n) n p (\[Beta] + p \[Eta] + \[Beta] \[Eta]))/(
kd (\[Beta] + p \[Eta]))]
BesselK[
2/(-3 + n), -((
rw^(3/2 - n/2) Sqrt[(
2^(1 + n) n p (\[Beta] + p \[Eta] + \[Beta] \[Eta]))/(
kd (\[Beta] + p \[Eta]))])/(-3 + n))] +
rw^2 Sqrt[(
2^(1 + n) n p (\[Beta] + p \[Eta] + \[Beta] \[Eta]))/(
kd (\[Beta] + p \[Eta]))]
BesselK[(
2 (-2 + n))/(-3 + n), -((
rw^(3/2 - n/2) Sqrt[(
2^(1 + n) n p (\[Beta] + p \[Eta] + \[Beta] \[Eta]))/(
kd (\[Beta] + p \[Eta]))])/(-3 + n))] +
2 (-1 + n) rw^((1 + n)/2)
BesselK[(-1 + n)/(-3 + n), -((
rw^(3/2 - n/2) Sqrt[(
2^(1 + n) n p (\[Beta] + p \[Eta] + \[Beta] \[Eta]))/(
kd (\[Beta] + p \[Eta]))])/(-3 + n))])), p, t, 16]

• You can read this post to learn how to format the code: mathematica.meta.stackexchange.com/q/1584/1871 Also, your code involves more than one simple mistake, please double check it. – xzczd Sep 16 '17 at 14:27
• The containing code has been edited. It looks better right now. I found the boundary condition is wrong. I change the boundary condition from Dirichlet to Neumann type. – LingLong Sep 17 '17 at 0:46
• What is rw` here? – zhk Sep 17 '17 at 1:48
• sorry for the typo. I revised it as R – LingLong Sep 17 '17 at 6:06
• I found that it can evaluate the numerical result when n less than 1. But n larger than 1 the MMA would not give the solution. – LingLong Sep 18 '17 at 12:01