PDE involving derivative at boundary, with a boundary condition at infinity

Question

I am trying to find the function $T(z,t)$ which solves this differential equation: $$\frac{\partial T}{\partial t}=\frac{\partial^2 T}{\partial z^2}+St\left[ \exp\left [ -\frac{\left( x_f-Ut\right )^2}{2\sigma^2} \right ] +\frac{\partial T}{\partial z}(0,t)\right]\frac{\partial T}{\partial z}$$ with the initial condition $$T(z,0)=T_0(z)$$ where $$T_0=T_c \exp (-a_1 z_c z)$$ and boundary conditions: $$T(0,t)=1$$ and $$T(\infty,t)=0$$ In the differential equation $$\frac{\partial T}{\partial z}(0,t)$$ is the derivative of the unknown function $T(z,t)$ calculated at $z=0$, while $x_f$ is calculated as the solution of the folowing equation: $$\exp\left ( -\frac{ x_f^2}{2\sigma^2} \right ) =-\frac{d T_0}{d z}(z=0)$$ The values of the constants are: $a_1=50$, $St=2$, $U=0.8$, $\sigma=0.4$, $Tc=1500$ and $z_c=0.001$. I tried the following code, even if I do not know if I have inserted correctly the term $\frac{\partial T}{\partial z}(0,t)$:

Tc = 1500
zc = 0.001
tc = 0.1
St = 2
U = 0.8
\[Sigma] = 0.4
a1=50;
T0[z_] := Exp[(-a1)*zc*z]
xf = NSolve[
   Exp[-(x^2/(2*\[Sigma]^2))] == -(D[T0[z], z] /. z -> 0) && x > 0, 
   x][[1, 1, 2]]
solu = NDSolve[{D[T[z, t], {t}] == 
    St*D[T[z, t], {z}]*(D[T[0, t], {z}] + 
        E^(-((xf - t*U)^2/(2*\[Sigma]^2)))) + 
           D[T[z, t], {z, 2}], T[0, t] == 1, T[190, t] == 0, 
   T[z, 0] == T0[z]}, T, {z, 0, 190}, {t, 0, 1000}]

I use $190$ to approximate $\infty$ here. Unfortunately I get an error of the type:

General::munfl: exp(-1.99511*10^6) is too small to represent as a normalized machine number; precision may be lost.

However the result seems to work and

D[T0[z], z] /. z -> 0

and

D[Evaluate[T[z, t] /. First[solu]], z] /. {z -> 0, t -> 0}

are quite similar, as they should. Indeed

(D[T0[z], z] /. 
   z -> 0) - (D[Evaluate[T[z, t] /. First[solu]], z] /. {z -> 0, 
    t -> 0})

is equal to $-0.0000357858$. But if I change $190$ to $700$ the previous difference becomes $-0.0036338$. This means that the result is very sensitive to the value used to approximate $\infty$. Is there any way to insert the condition at $\infty$ in NDSolve without using a numerical approximation, also because I am interested in studying how $T(z,t)$ varies with $a_1$, and I noticed that varying $a_1$ I need to vary the numerical approximation of $\infty$ to have good results. Many thanks. (Could you please check if I wrote the Mathematica code for the differential equation correctly).

Answer 1

Perhaps NestList, assuming D[T[0, t], {z}] is known from pervious iteration, helps to find a solution.

sol = NestList[
NDSolveValue[{D[T[z, t], {t}] ==St*D[T[z, t], {z}]*(
Derivative[1, 0][#] [0, t] + 
E^(-((xf - t*U)^2/(2*\[Sigma]^2)))) + D[T[z, t], {z, 2}], 
T[0, t] == 1, T[190, t] == 0, T[z, 0] == T0[z]}, 
T, {z, 0, 190}, {t, 0, 1000}] &, 0 &, 10]

This gives a solution which seems to converge:

Row[{Plot3D[{ sol[[-1]][z, t]} , {z, 0, 190}, {t, 0, 1000},ImageSize -> 200],
Plot3D[{sol[[-2]][z, t] - sol[[-1]][z, t]} , {z, 0, 190}, {t, 0,1000}, ImageSize-> 200]}]

Part D[T[0, t], {z}] converges too:

Plot[Map[Derivative[1, 0][Function[{z, t}, #[z, t]]][0, t] &,sol[[-5 ;;]]], {t, 0, 1000}, Evaluated -> True]

By the way: You might omit the infinity bc T[190, t] == 0 and get the same result!

Hope it helps!

Answer 2

Let me extend my comment to an answer.

First of all, D[T[0, t], {z}] is obviouly wrong. Since there's literally no z in T[0, t], D[T[0, t], {z}] will evaluate to 0:

"But NDSolve is so powerful! There may be some magic in NDSolve so D[T[0, t], {z}] will be parsed as I expect!" <- If this is what you're thinking, please notice NDSolve doesn't have HoldAll/HoldFirst/… attribute so its arguments are all evaluated before passed into NDSolve i.e. NDSolve doesn't even have a chance to analyze the D[T[0, t], {z}] term. (A more detailed analysis can be found here. )

"OK, then how to circumvent?" Currently NDSolve cannot handle this problem directly. Picard iteration as shown in Ulrich's answer is a way to go. Discretizing the PDE in $z$ direction with pdetoode is another:

domain = {lb, rb} = {0, 190};
tend = 1000;
With[{T = T[z, t]}, 
    {eq, ic, bc} = 
     {D[T, t] == St D[T, z] (term + E^(-((xf - t U)^2/(2 σ^2)))) + D[T, {z, 2}], 
      T == T0[z] /. t -> 0, 
      {T == 1 /. z -> lb, T == 0 /. z -> rb}};
  termvalue = D[T, z] /. z -> lb];

points = 25; difforder = 4; grid = Array[# &, points, domain];
(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)
ptoofunc = pdetoode[T[z, t], t, grid, difforder];
del = #[[2 ;; -2]] &;

ode = del@ptoofunc@eq /. term -> ptoofunc@termvalue;
odeic = ptoofunc@ic // del;
odebc = ptoofunc@bc;

sollst = NDSolveValue[{ode, odeic, odebc}, T /@ grid, {t, 0, tend}]; // AbsoluteTiming
sol = rebuild[sollst, grid, 2];

Plot3D[sol[z, t], {z, lb, rb}, {t, 0, tend}]

Plot[termvalue /. T -> sol // Evaluate, {t, 0, tend}]