Solving the heat equation on the semi-infinite rod

Question

Cross posted in scicomp.SE.

---

I want to test the solution which is given below is right by Mathematica.

Please look the post in mathstackexhange

or

Please look below.

Question: Solve the following heat equation on the semi-infinite rod

$\qquad u_t=ku_{xx}$

where $x,t>0$ and

$\qquad u_x(0,t) =0$ and $u(x,0)=\begin{cases} 1, & 0 < x <2 \\ 0, & 2\leq x \end{cases}$

with proper Fourier transform.

Answer

We found the following answer

$ \qquad u(x,t) = \frac{2}{\pi}\int_{0}^{\infty}e^{-s^2 t}\frac{1-\cos(2s)}{s}\cos(sx)ds. $

Code

But I am not sure the solution is right. I am not capable of testing it in Mathematica.

Could you help me?

Answer 1

Personally I think the problem is interesting, so let me extend my comments to an answer. First of all, DSolve can solve OP's problem straightforwardly (in Mathematica 10.3 or higher, if I remember correctly):

With[{u = u[t, x]}, 
 eq = D[u, t] == k D[u, x, x];
 ic = u == Piecewise[{{1, 0 < x < 2}}] /. t -> 0;
 bc = D[u, x] == 0 /. x -> 0;]

asol = DSolveValue[{eq, ic, bc}, u, {t, x}, Assumptions -> {x > 0, k > 0}];
asol[t, x]
(* 1/2 (-Erf[(-2 + x)/(2 Sqrt[k] Sqrt[t])] + Erf[(2 + x)/(2 Sqrt[k] Sqrt[t])]) *)

Remark

There seems to be a bug in DSolve in v11.2.0.

DSolve[{eq, ic, bc}, u[t, x], {t, x}]

will return unevaluated.

As one can see, DSolve expresses the solution with Erf, so it's not immediately clear whether OP's solution is correct or not, and Mathematica's functions for simplifying also doesn't work well in this case, so let's obtain the analytic solution with another approach, that is, making use of Fourier cosine transform to eliminate the derivative of $x$:

fct = FourierCosTransform[#, x, s] &;

tset = Map[fct, {eq, ic}, {2}] /. Rule @@ bc /. 
  HoldPattern@FourierCosTransform[a_, __] :> a

tsol = u[t, x] /. DSolve[tset, u[t, x], t][[1]]
(* (E^(-k s^2 t) Sqrt[2/π] Sin[2 s])/s *)

Remark

I've made the transform on the PDE in a quick way, for a more general approach, check this post.

InverseFourierCosTransform has difficulty in transforming tsol, but it doesn't matter because the integral form is just what we want. By checking the formula of inverse Fourier cosine transform, we find the solution should be

$$u(t,x)=\sqrt{\frac{2}{\pi }} \int_0^{\infty } \frac{e^{-k s^2 t} \sqrt{\frac{2}{\pi }} \cos (s x) \sin (2 s)}{s} \, ds$$

It's apparently different from the one in your question, and numeric calculation shows this solution is the same as the one given by DSolve, so the one in your question is wrong.

Finally, a illustration for the solution:

Plot3D[asol[t, x] /. k -> 1 // Evaluate, {x, 0, 4}, {t, 0, 10}]

---

Update

Inspired by Ars3nous' comment below, I noticed InverseFourierCosTransform can actually transform tsol. We just need a proper assumption:

InverseFourierCosTransform[tsol, s, x, Assumptions -> k > 0]
(* 1/2 (-Erf[(-2 + x)/(2 Sqrt[k t])] + Erf[(2 + x)/(2 Sqrt[k t])]) *)

Apparently it's the same as asol.