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
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


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).

  • 3
    $\begingroup$ D[T[0, t], {z}] is obviouly wrong. (Just execute it separately and observe! ) NDSolve cannot handle this problem directly, you need to discretize yourself. $\endgroup$
    – xzczd
    Jan 1, 2023 at 1:04

2 Answers 2


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]}]

enter image description here

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]

enter image description here

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

Hope it helps!

  • $\begingroup$ Thank you agian. I was surprised by discovering that NDSolve works also without the bc at infinity, since the problem should be ill-posed in this way. Indeed, the kernel issues a warning in this case: an insufficient number of boundary conditions have been specified for the direction of independent variable z. However, the solution produced without the abovementioned bc is practically identical to the one you get with the bc. Do Mathematica automatically add a bc to close the problem? Is this the meaning of the following warning: artificial boundary effects may be present in the solution. $\endgroup$
    – umby
    Jan 4, 2023 at 10:35
  • 1
    $\begingroup$ @umby No, NDSolve doesn't add any b.c., it simply uses one-sided difference formula to close the system. Please notice the bcart warning is severe and you should avoid it. For more info, read: mathematica.stackexchange.com/q/73961/1871 $\endgroup$
    – xzczd
    Jan 4, 2023 at 10:59

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:

enter image description here

"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}]

enter image description here

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

enter image description here

  • 2
    $\begingroup$ To the downvoter, I am interested in what's wrong with my answer, would you please elaborate. I'm not trying to complain here, I just want improve my answer if possible. $\endgroup$
    – xzczd
    Jan 4, 2023 at 1:12
  • $\begingroup$ Dear xzczd, dear Ulrich Neumann, many thanks! your answers really solved my problem, I'm using them both and having fun in solving my pde also changing the initial condition. For this reason I accept them both as a good solution for my problem. $\endgroup$
    – umby
    Jan 4, 2023 at 10:20
  • $\begingroup$ Dear @xzczd , just to be sure, what is difforder? why difforder is 4 in this case? $\endgroup$
    – umby
    Jul 19, 2023 at 10:12
  • 1
    $\begingroup$ @user58367 It essentially sets the value of DifferenceOrder option of NDSolve`FiniteDifferenceDerivative, on which pdetoode is built. DifferenceOrder -> 4 is a common choice, another common choice is 2. Some discussion for this option can be found in the built-in tutorial The Numerical Method of Lines. $\endgroup$
    – xzczd
    Jul 19, 2023 at 10:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.