2
$\begingroup$

I am wondering what is the best way to solve numerically a initial-boundary value problem of heat equation like this (surely I've made a lot of mistakes in the code):

NDSolve[{4 D[T[x, t], t] == D[T[x, t], x, x] , T[x, -Infinity] == 0, T[Infinity, t] == 0, -Derivative[1, 0][T][0, t] == 2 Sqrt[Pi] Exp[-(t^2/2)] }, T, {x, 0, 5}, {t, 0, 5}]

As we can see, there's one i.c. at $-\infty$:

T[x, -Infinity] == 0

one b.c. at $x=\infty$:

T[Infinity, t] == 0

and one b.c. as the derivative of the unknown with respect to $x$ at $x=0$ equal to the function $\exp(-t^2/2)$ (i.e. $-\left.\frac{\partial T}{\partial x}\right|_{x=0}=2\sqrt{\pi} \exp(-t^2/2)$)

-Derivative[1, 0][T][0, t] == 2 Sqrt[Pi] Exp[-(t^2/2)]

The analytical solution is:

T[x_, t_] := NIntegrate[ 1/Sqrt[τ] Exp[-((t - τ)^2/2)] Exp[-(x^2/τ)], {τ, 0, Infinity}]

Just to check the solution:

With[{trange = {-2, 0.5, 3}}, Show[{ListPlot[ Table[{x, NIntegrate[ 1/Sqrt[τ] Exp[-((t - τ)^2/ 2)] Exp[-(x^2/τ)], {τ, 0, Infinity}]}, {t, trange}, {x, 0, 4, 0.05}], PlotLegends -> {"t = -2", "t = 0.5", "t = 3"}], Table[Plot[ NIntegrate[ 1/Sqrt[τ] Exp[-((t - τ)^2/ 2)] Exp[-(0^2/τ)], {τ, 0, Infinity}] - 2 Sqrt[Pi] Exp[-(t^2/2)] x, {x, 0, 0.7}], {t, trange}]}, Frame -> True, PlotRange -> All, FrameLabel -> {x, "T[x,t]"}]]

enter image description here

Thanks in advance.

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ @umby Would you please explain the different boundary conditions, especially T[x, -Infinity] == 0 . Thanks. $\endgroup$ Commented Oct 9 at 17:39
  • 1
    $\begingroup$ Can you add a bit more background info for the problem? That'll make the problem more interesting, this type of problem is never discussed in this site AFAIK. $\endgroup$ Commented Oct 10 at 2:08
  • $\begingroup$ @UlrichNeumann, I think T[x, -Infinity] == 0 can be interpreted as an initial condition. The BCs are T[Infinity,t] == 0 and the condition on the derivative: -Derivative[1, 0][T][0, t] == 2 Sqrt[Pi] Exp[-(t^2/2)]. Thanks. $\endgroup$ Commented Oct 10 at 7:54
  • $\begingroup$ @xzczd, basically it is a semi-infinite body heated by a Gaussian source at its surface in 1D. $\endgroup$ Commented Oct 10 at 8:01
  • 1
    $\begingroup$ 10 correct digits: NIntegrate[(Exp[(-(1/2)) (t - \[Tau])^2] Exp[-(x^2/\[Tau])])/ Sqrt[\[Tau]] /. t -> 1 /. x -> 1, {\[Tau], 0, Infinity}] - Sum[((-1)^n 2^(-1 - n) E^(-(t^2/2)) t^m x^( 4 n) ((2^(1/4 + m/2) Gamma[3/4 - m/2] Gamma[1/4 + m/2])/( Gamma[3/4 - m/2 + n] Gamma[1 + 2 n]) + x ((2^(-(1/4) + m/2) x Gamma[1/4 - m/2] Gamma[3/4 + m/2])/( Gamma[5/4 - m/2 + n] Gamma[2 + 2 n]) - ( 2 x^(2 m) Gamma[1/2 - m] Gamma[1/2 + m])/( Gamma[1 + n] Gamma[3/2 + m + 2 n]))))/Gamma[1 + m] /. t -> 1 /. x -> 1, {n, 0, 20}, {m, 0, 20}] $\endgroup$ Commented Oct 10 at 10:03

3 Answers 3

Reset to default
5
$\begingroup$

In many cases the solution of heat equation decays fast in time and space, so using a large enough number to approximate infinity often works well, and luckily your problem seems to be the case:

tend = 10.; inf = 5.; nsol = NDSolveValue[{4 D[T[x, t], t] == D[T[x, t], x, x], T[x, -inf] == 0, T[inf, t] == 0, -Derivative[1, 0][T][0, t] == 2 Sqrt[Pi] Exp[-(t^2/2)]}, T, {x, 0, inf}, {t, -inf, tend}]; (* Alternatively: *) (* nsol = NDSolveValue[{4 D[T[x, t], t] == D[T[x, t], x, x] + NeumannValue[2 Sqrt[Pi] Exp[-(t^2/2)], x == 0], T[x, -inf] == 0, T[inf, t] == 0}, T, {x, 0, inf}, {t, -inf, tend}]; *) Tsol[x_, t_] := NIntegrate[ (Exp[(-(1/2)) (t - τ)^2] Exp[-(x^2/τ)])/ Sqrt[τ], {τ, 0, Infinity}] lst = Table[{x, t, Tsol[x, t]}, {x, 0., inf, 1/2}, {t, -inf, tend, 1/2}]; ListPointPlot3D@lst~Show~ Plot3D[nsol[x, t], {x, 0, inf}, {t, -inf, tend}, PlotRange -> All]

enter image description here

Remark

  1. The problem would be much harder if you asked how to deduce the symbolic solution with the help of Mathematica.

  2. It's worth pointing out that, when verifying the PDE or b.c. at $x=0$, one cannot directly substitute in the analytic solution, because in this case integration and differentiation is not interchangable, but D or Derivative goes into the integral automatically. Quick example:

    Tsol2[x_, t_] := Inactive[Integrate][(Exp[-(1/2) (t - τ)^2] Exp[-(x^2/τ)])/Sqrt[τ], {τ, 0, ∞}] D[Tsol2[x, t], x] /. x -> 0 (* Inactive[Integrate][0, {τ, 0, ∞}] <- Obviously wrong. *)

    To obtain the correct result, one possible way is to use ND:

    Clear[Tsol3]; Tsol3[x_?NumericQ, t_] := NIntegrate[(Exp[-(1/2) (t - τ)^2] Exp[-(x^2/τ)])/Sqrt[τ], {τ, 0, ∞}] << NumericalCalculus` ND[Tsol3[x, 1], x, 0] (* -2.1501 *)
Improve this answer
$\endgroup$
7
  • 1
    $\begingroup$ @xzczd Nice answer. I only checked the residuum of the pde for the proposed analytical answer,that deviates significantly from zero. Don't know why $\endgroup$ Commented Oct 10 at 8:09
  • 1
    $\begingroup$ @umby Nothing special here, just repeat the analysis I've done in e.g. mathematica.stackexchange.com/a/245309/1871 , you'll know the correct Neumann value. $\endgroup$ Commented Oct 10 at 8:21
  • 1
    $\begingroup$ @Ulrich Just noticed there's a typo in your code, Derivative[0, 2][T][x, t] should be Derivative[0, 1][T][x, t]. After correcting this i.e. using resi = Function[{x, t}, 4 Derivative[0, 1][Tsol][x, t] - Derivative[2, 0][Tsol][x, t]], the result looks fine. $\endgroup$ Commented Oct 10 at 8:26
  • 1
    $\begingroup$ @xzczd Upss, thanks! $\endgroup$ Commented Oct 10 at 12:25
  • 2
    $\begingroup$ @xzczd The finite element ( solT = NDSolveValue[{4 D[T[x, t], t] == D[T[x, t], x, x] + NeumannValue[2 Sqrt[Pi] Exp[-(t^2/2)], x == 0], T[x, -inf] == 0, T[inf, t] == 0}, T, {x, 0, inf}, {t, -inf, inf}, Method -> {"MethodOfLines", "TemporalVariable" -> t, "SpatialDiscretization" -> {"FiniteElement"}}] ) solution also fits the region x==0 very well. Only the point x==0,t==0 shows residual unequal to zero, which can be explained by the inconsistent initial conditions. $\endgroup$ Commented Oct 12 at 9:58
2
$\begingroup$

modified

To long for a comment:

Your "analytical" solution

T = Function[{x, t},NIntegrate[1/Sqrt[\[Tau]] Exp[-((t - \[Tau])^2/2)] Exp[-(x^2/\[Tau])], {\[Tau], 0, Infinity},Method -> "LocalAdaptive"]]

should satisfy the PDE to a good approximation

resi = Function[{x, t},4 Derivative[0, 1][T][x, t] - Derivative[2, 0 ][T][x, t]]; Plot3D[resi[x, t], {x, 0, 1}, {t, 0, 1}, PlotRange -> All,AxesLabel -> {x, t, residual}]

enter image description here

But plot now shows a residual which equals zero, except for x==0.

Improve this answer
$\endgroup$
8
  • $\begingroup$ Ulrich, I’m afraid I forgot a factor of 4; I’ve corrected the post, and I hope it’s fine now. Please accept my sincerest apologies. $\endgroup$ Commented Oct 9 at 16:47
  • $\begingroup$ It should be: T = Function[{x, t}, NIntegrate[ 1/Sqrt[\[Tau]] Exp[-((t - \[Tau])^2/ 2)] Exp[-(x^2/\[Tau])], {\[Tau], 0, Infinity}, Method -> "LocalAdaptive"]] resi = Function[{x, t}, 4 Derivative[0, 1][T][x, t] - Derivative[2, 0][T][x, t]]; Plot3D[resi[x, t], {x, 0, 1}, {t, 0, 1}], it's the first derivative with restect to t. $\endgroup$ Commented Oct 9 at 16:54
  • 1
    $\begingroup$ @umby But the pde does not fulfill this function either! $\endgroup$ Commented Oct 9 at 17:37
  • 1
    $\begingroup$ @UlrichNeumann This time it's right, see my answer for more info. $\endgroup$ Commented Oct 10 at 2:07
  • 1
    $\begingroup$ @xzczd Thanks for your correction, I modified my answer. Now the residual vanishs, except near x=0 . Any idea why? $\endgroup$ Commented Oct 10 at 12:35
1
$\begingroup$

I still don't understand the constraints T[x, -Infinity] == 0,T[ Infinity,t] == 0

Without these constraints we can find an numerical solution assuming NeumannValue[0 ,x==inf] and NeumannValue[2 Sqrt[Pi] Exp[-(t^2/2)], x == 0] 

inf = 20; (* nuemrical infinity*) tmax = 20; z =NDSolveValue[{4 Derivative[0, 1][T][x, t] -Derivative[2, 0][T][x, t] == NeumannValue[2 Sqrt[Pi] Exp[-(t^2/2)], x == 0] , T[x, 0] == 0 }, T, {x, 0, inf}, {t, 0, tmax}, Method -> {"MethodOfLines", "TemporalVariable" -> t, "SpatialDiscretization" -> {"FiniteElement" }}] Plot3D[z[x, t], {x, 0 m, inf}, {t, 0, tmax}, PlotRange -> All,MeshFunctions -> (#2 &)]

enter image description here

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.