# Solving the heat equation on the semi-infinite rod

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

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.

$\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?

• I don't we can do much to help without having the code you used to get the answer. If you didn't use Mathematica to solve your problem, then this question is inappropriate -- I would mean you are asking us to both write the code and verify the solution for you. Oct 17 '17 at 17:34
• Well, solving this in Mathematica is quite straightforward, just check the document of DSolve. Anyway, the solution is 1/2 (-Erf[(-2 + x)/(2 Sqrt[k t])] + Erf[(2 + x)/(2 Sqrt[k t])]). A quick test shows the solution in your question seems to be wrong. Oct 18 '17 at 4:44
• Further check shows that, if one wants to express the solution as integration, then it should be $\int_0^{\infty } \frac{\sqrt{\frac{2}{\pi }} \sin (2 w) e^{-k t w^2} \cos (w x)}{w} \, dw$ Oct 18 '17 at 12:08
• Guys, personally I suggest not to close this post, though it's a… "give me the code" question, the problem is interesting, I think. Oct 18 '17 at 12:13
• Please do not cross-post within SE sites. Choose one site and delete the questions from the others. meta.stackexchange.com/q/64068/164803 Oct 19 '17 at 12:17

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.

• Following @xzczd amazing answer using cosine transform, with k=1 (you can always scale out k in time) In Mathematica 9.0, I get using InverseFourierCosineTransform, $u(x,t)=\frac{1}{2} \left(-2 x \text{erfc}\left(\frac{x}{2 \sqrt{t}}\right)+\text{erfc}\left(\frac{x+2}{2 \sqrt{t}}\right)+\text{erfc}\left(-\frac{x-2}{2 \sqrt{t}}\right)+\frac{4 \sqrt{t} e^{-\frac{x^2}{4 t}}}{\sqrt{\pi }}-2\right)$ which has two different terms along the mentioned answer. But this solution also does not satisfy boundary conditions. I wonder what is the discrepency for. Nov 5 '17 at 12:18
• @Ars3nous Are you in v9.0.0 or v9.0.1? I just tested in v9.0.1, Win10 64bit, with InverseFourierCosTransform[tsol /. k -> 1, s, x] I got 1/2 (Erf[(2 - x)/(2 Sqrt[t])] + Erf[(2 + x)/(2 Sqrt[t])]), which is consistent with asol. Nov 5 '17 at 12:38