Heat equation with nonlinear boundary condition involving time-derivative

The governing equation is shown as follows: I first try to employ the NDSolve, but it seems that Mathematica can not handle the fourth boundary condition.

Therefore, I rewrite the code with finite difference method, but still have problem

Mr = 100; Mz = 10; Lr = 10; Lz = 1;
γ = 0.01;
α = 1;
κ = 1;
dr = Lr/Mr;
dz = Lz/dz;
dt = 0.1;
r[i_] := i*dr;
z[j_] := j*dz;

(*i.c.*)
u[i_, j_, 0] := 0;
(*b.c.*)
u[1, j_, k_] := u[2, j, k] + 2;
u[Mr, j_, k_] := u[Mr - 1, j, k];
u[i_, 1, k_] := u[i, 2, k];
u[i_, Mz, k_] :=
If[k == 1, 0,
1/dz^2 (dt dz γ u[i, -1 + Mz, -1 + k] +
dt α γ u[i, -1 + Mz, -1 + k]^2 +
dz^2 u[i, Mz, -1 + k] - dt dz γ u[i, Mz, -1 + k] -
2 dt α γ u[i, -1 + Mz, -1 + k] u[i, Mz, -1 + k] +
dt α γ u[i, Mz, -1 + k]^2)];

u[i_, j_, k_] :=
u[i, j, k] =
1/(dr^2 dz^2 i) (dt dz^2 i u[-1 + i, j, -1 + k] +
dr^2 dt i u[i, -1 + j, -1 + k] - dt dz^2 u[i, j, -1 + k] -
2 dr^2 dt i u[i, j, -1 + k] + dr^2 dz^2 i u[i, j, -1 + k] -
2 dt dz^2 i u[i, j, -1 + k] + dr^2 dt i u[i, 1 + j, -1 + k] +
dt dz^2 u[1 + i, j, -1 + k] + dt dz^2 i u[1 + i, j, -1 + k]);

u[1, Mz, 1]

However, I receive the following error message

$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>

$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >> General::stop: Further output of$RecursionLimit::reclim will be suppressed during this calculation. >>

The problem may be induced by the imposion of nonlinear and time-derivative and nonlinear boundary condition. Any idea about how to solve my problem? This boundary condition has bothered me for long time.

I re code the provided code as follows

Clear[fdd, pdetoode, tooderule]
fdd[{}, grid_, value_, order_] := value;
fdd[a__] := NDSolveFiniteDifferenceDerivative@a;

pdetoode[funcvalue_List, rest__] :=
pdetoode[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[],
rest];
pdetoode[{func__}[var__], rest__] :=
pdetoode[Alternatives[func][var], rest];
pdetoode[rest__, grid_?VectorQ, o_Integer] :=
pdetoode[rest, {grid}, o];

pdetoode[func_[var__], time_, {grid : {__} ..}, o_Integer] :=
With[{pos = Position[{var}, time][[1, 1]]},
With[{bound = #[[{1, -1}]] & /@ {grid},
pat = Repeated[_, {pos - 1}],
spacevar = Alternatives @@ Delete[{var}, pos]},
With[{coordtoindex =
Function[coord,
Piecewise[{{1, # === #2[]}, {-1, # === #2[[-1]]}},
All] &, {coord, bound}]]},
tooderule@
Flatten@{((u : func) |
Derivative[dx1 : pat, dt_, dx2___][(u : func)])[x1 : pat,
t_, x2___] :> (Sow@coordtoindex@{x1, x2};

fdd[{dx1, dx2}, {grid},
Outer[Derivative[dt][u@##]@t &, grid],
"DifferenceOrder" -> o]),
inde : spacevar :>
With[{i = Position[spacevar, inde][[1, 1]]},
Outer[Slot@i &, grid]]}]]];

tooderule[rule_][pde_List] := tooderule[rule] /@ pde;
tooderule[rule_]@Equal[a_, b_] :=
Equal[tooderule[rule][a - b], 0] //.
eqn : HoldPattern@Equal[_, _] :> Thread@eqn;
tooderule[rule_][expr_] := #[[Sequence @@ #2[[1, 1]]]] & @@
Reap[expr /. rule]

Clear@pdetoae;
pdetoae[funcvalue_List, rest__] :=
pdetoae[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[], rest];
pdetoae[{func__}[var__], rest__] :=
pdetoae[Alternatives[func][var], rest];

pdetoae[func_[var__], rest__] :=
Module[{t},
Function[
pde, #[pde /. {Derivative[d__][u : func][inde__] :>
Derivative[d, 0][u][inde, t], (u : func)[inde__] :>
u[inde, t]}] /. (u : func)[i__][t] :> u[i]] &@
pdetoode[func[var, t], t, rest]];
\[Gamma] = 100;
\[Alpha] = 0;
\[Kappa] = 1;
R = 10;
Z = 1;
eps = 10^-1;
tend = 1000;

eq = With[{u = u[r, Sqrt[\[Kappa]] z, t]},
Laplacian[u, {r, th, z}, "Cylindrical"] == D[u, t] /.
Sqrt[\[Kappa]] z -> z];

ic = u[r, z, 0] == 0;
bc = With[{u = u[r, z, t]}, {D[u, r] == -2/r /. r -> eps,
u == 0 /. r -> R, D[u, z] == 0 /. z -> 0, D[u, z] == 0 /. z -> Z}]

domain@r = {eps, R};
domain@z = {0, Z};
points@r = 50;
points@z = 50;
difforder = 2;
(grid@# = Array[# &, points@#, domain@#]) & /@ {r, z};

(*Definition of pdetoode isn't included in this post,please find it \
in the link above.*)
ptoofunc = pdetoode[u[r, z, t], t, grid /@ {r, z}, difforder];

delbothside = #[[2 ;; -2]] &;

ode = delbothside /@ delbothside@ptoofunc@eq;

odeic = delbothside /@ delbothside@ptoofunc@ic;
odebc = MapAt[delbothside, ptoofunc@bc, {{1}, {2}}];

sollst = NDSolveValue[{ode, odeic, odebc},
Outer[u, grid@r, grid@z] // Flatten, {t, 0, tend}];

sol = ListInterpolation[
Partition[DeveloperToPackedArray@#["ValuesOnGrid"] & /@ #,
points@z], {grid@r, grid@z, #[]["Coordinates"][]}] &@
sollst; // AbsoluteTiming

Manipulate[
Plot3D[sol[r, z, t], {r, eps, R}, {z, 0, Z},
PlotRange -> {-1, 10}], {t, 0, tend}]

The equation indicate that the upper and under boundary conditions are subjected to no-flow condition. And it can easily be solved and obtain an analytical solution.

To compare the analytical solution (or semi-analytical), the solution can be expressed as

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/t Sum[ Vi[n, i]*F /. s -> i Log/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/t Sum[
Vi[n, i]*F /. s -> i Log/t , {i, 1, n} ] ]]]  // N;
s0[r_, z_, t_] :=
NIntegrate[
Re[Stehfest[(
2 E^(-((Sqrt[a^2 + p] z)/
Sqrt[\[Kappa]])) (-E^((Sqrt[a^2 + p]/Sqrt[\[Kappa]]))
p \[Gamma] -
E^((Sqrt[a^2 + p] (1 + 2 z))/Sqrt[\[Kappa]]) p \[Gamma] +
E^((Sqrt[a^2 + p] z)/
Sqrt[\[Kappa]]) (p \[Gamma] - Sqrt[a^2 + p] Sqrt[\[Kappa]]) +
E^((Sqrt[a^2 + p] (2 + z))/
Sqrt[\[Kappa]]) (p \[Gamma] +
Sqrt[a^2 + p] Sqrt[\[Kappa]])))/(
p (a^2 +
p) ((1 + E^((2 Sqrt[a^2 + p])/
Sqrt[\[Kappa]])) p \[Gamma] + (-1 + E^((2 Sqrt[a^2 + p])/
Sqrt[\[Kappa]])) Sqrt[a^2 + p] Sqrt[\[Kappa]])), p, t, 6]]*
BesselJ[0, a (r)]*a, {a, 0, \[Infinity]},
Method -> {"LevinRule", "LevinFunctions" -> {BesselJ}},
MaxRecursion -> 40];

with \[Gamma] = 0; I however compare the numerical results and analytical one. It seems that the the numerical solution obtain a wrong result. The code and the plot of comparison are shown as follows

\[Gamma] = 0; TT1 = Table[{t, s0[0.9, 0.5, t]}, {t, 0.1, 1, 0.1}];
TT2 = Table[{t, s0[0.9, 0.5, t]}, {t, 1, 10, 1}];
TT3 = Table[{t, s0[0.9, 0.5, t]}, {t, 10, 100, 10}];
TT4 = Table[{t, s0[0.9, 0.5, t]}, {t, 100, 1000, 100}];
C0 = Join[TT1, TT2, TT3, TT4];
TT1 = Table[{t, sol[1, 0.5, t]}, {t, 0.1, 1, 0.1}];
TT2 = Table[{t, sol[1, 0.5, t]}, {t, 1, 10, 1}];
TT3 = Table[{t, sol[1, 0.5, t]}, {t, 10, 100, 10}];
TT4 = Table[{t, sol[1, 0.5, t]}, {t, 100, 1000, 100}];
Cn = Join[TT1, TT2, TT3, TT4]; The figure is u versus t and blue and yellow lines represent the value predicted by numerical and analytical solutions, respectively. I an trying to figure out what happen to the numerical method.

• As mentioned in the comment below, please add a complete example for producing the graph. – xzczd May 2 '17 at 10:07
• @xzczd The full code is given. – LingLong May 2 '17 at 10:32
• The γ in 2 approaches have different definitions? – xzczd May 2 '17 at 11:17
• The same meaning. But I need to let it become zero to get the same solution form to compare to numerical result. Because the analytical solution is derived from the use of the nonlinear boundary condition and let gamma zero to reduce the upper boundary condition (i.e., the nonlinear boundary) to no-flow one. – LingLong May 2 '17 at 11:37
• Have you read this (BTW here's the unfinished Chinese edition of it) and this book? – xzczd May 2 '17 at 14:32

NDSolve has trouble in handling the last b.c., so let's help it a bit by discretizing the PDE and corresponding i.c. and b.c. to a set of DAE.

First, interpreting your equation set to Mathematica code:

γ = 1/100;
α = 1;
κ = 1;
R = 10;
Z = 1;
eps = 10^-1;
tend = 1;

eq = With[{u = u[r, Sqrt[κ] z, t]},
Laplacian[u, {r, th, z}, "Cylindrical"] == D[u, t] /. Sqrt[κ] z -> z];

ic = u[r, z, 0] == 0;
bc = With[{u = u[r, z, t]},
{D[u, r] == -2/r /. r -> eps,
u == 0 /. r -> R,
D[u, z] == 0 /. z -> 0,
-1/γ eq[] == D[u, z] - α D[u, z]^2 /. z -> Z}]

Notice I've substitute the PDE to the last b.c., because it turns out that, NDSolve can't handle the original one very well even after discretization.

Then it's time for discretization, I'll use pdetoode for the task:

domain@r = {eps, R};
domain@z = {0, Z};
points@r = 25;
points@z = 25;
difforder = 4;
(grid@# = Array[# &, points@#, domain@#]) & /@ {r, z};

(* Definition of pdetoode isn't included in this post,
please find it in the link above. *)
ptoofunc = pdetoode[u[r, z, t], t, grid /@ {r, z}, difforder];

delbothside = #[[2 ;; -2]] &;

ode = delbothside /@ delbothside@ptoofunc@eq;

odeic = delbothside /@ delbothside@ptoofunc@ic;
odebc = MapAt[delbothside, ptoofunc@bc, {{1}, {2}}];

sollst = NDSolveValue[{ode, odeic, odebc}, Outer[u, grid@r, grid@z], {t, 0, tend}];

Finally, rebuild the solution from the discretized one and plot:

sol = rebuild[sollst, grid /@ {r, z}, -1]; // AbsoluteTiming

Manipulate[Plot3D[sol[r, z, t], {r, eps, R}, {z, 0, Z}, PlotRange -> {-1, 10}],
{t, 0, tend}] • What's the amazing work! thank you! I need time to understand your code. – LingLong Apr 30 '17 at 9:31
• (+1) BTW, I still didn't understand how you input that weird bc. It's (MMM) now zhk. – zhk Apr 30 '17 at 9:49
• this bc represents the drainage from the top. Some one linearize the moving bc like this. However, i m trying to input the linearized one (i.e., alpha is equal to zero), the numerical results can not match the semi-analytical solution (the equations with alpha = 0 can be solved). It produce the negative values somewhere. – LingLong Apr 30 '17 at 10:00
• @zhk The left hand side (RHS) of the weird b.c. is $\partial{u}/\partial{t}$, which is the same as the RHS of the PDE, so I substitute the PDE to that b.c.. (Just evaluate eq[] separately and see what will happen. ) – xzczd Apr 30 '17 at 10:39
• @JOJO I believe you mean something like that mentioned in this post?: mathematica.stackexchange.com/q/10055/1871 – xzczd Apr 30 '17 at 10:41