Inhomogeneous Heat Equation Formula Not Satisfying IVP

Question

From Wikipedia,

However, when I plug the formula of $u(x,t)$ into Mathematica 12.1, it doesn't seem to satisfy the PDE (does not give $0$):

Here is the code:

G[x_, t_] := 1/Sqrt[4*Pi*k*t]*Exp[-x^2/(4*k*t)]

u[x_, t_] := Integrate[G[x - y, t - s]*f[y, s], {s, 0, t}, {y, -Infinity, Infinity}]

FullSimplify[D[u[x, t], t] - k*D[u[x, t], {x, 2}] - f[x, t], Assumptions -> t > 0 && k > 0]

Answer 1

The answer by Mathematica is correct

Clear["Global`*"]
pde = D[u[x, t], t] == k D[u[x, t], {x, 2}] + f[x, t];
ic = u[x, 0] == 0;
sol = DSolve[{pde, ic}, u[x, t], {x, t}];
sol /. {K[1] -> y, K[2] -> s}

$$ \left\{\left\{u(x,t)\to \frac{\int _0^t\int _{-\infty }^{\infty }\sqrt{\frac{k}{t-s}} f(y,s) e^{-\frac{(x-y)^2}{4 k (t-s)}}dyds}{2 \sqrt{\pi } k}\right\}\right\} $$

Which is the same as one you show

Verfication

The hint below given by xzczd made me think that if $f(x,t)$ was given, i.e. a specific function, then Mathematica should be able to verify the book solution.

And this is indeed the case. Tried few random $f(x,t) and Mathematica can verify it now, giving True. It does take few seconds to do, depending how complicated $f(x,t)$ is.

Clear["Global`*"];

moveconst is from how to simplify symbolic integration thanks to celtschk

moveconst[
   x_] := (x /. 
    Integrate[factor_ expr_, {var_, min_, max_}] /; 
      FreeQ[factor, var] :> factor Integrate[expr, {var, min, max}]);

pde = D[u[x, t], t] == k*D[u[x, t], {x, 2}] + f[x, t];
G[x_, t_] := 1/Sqrt[4*Pi*k*t]*Exp[-x^2/(4*k*t)];
bookAnswer = 
  Integrate[G[x - y, t]*g[y], {y, -Infinity, Infinity}] + 
   Integrate[
    G[x - y, t - s]*f[y, s], {s, 0, t}, {y, -Infinity, Infinity}];
f[x_, t_] := x*t;
sol = u -> Function[{x, t}, Evaluate@bookAnswer];
result = pde /. sol;
Simplify[result, TransformationFunctions -> {Automatic, moveconst}]

(*True*)

Tried f[x_, t_] := Sin[t]*x; and few others, all give True.

Answer 2

I'd like to extend my comment to an answer. As mentioned above, OP's attempt fails mainly because D calculates based on Leibniz integral rule but sadly this rule doesn't apply in this case, because the integrand is no longer continuous when $s=t$.

Still, we can verify the solution using Mathematica, with a little manual analysis. First of all, we introduce a positive $\epsilon$ to the solution:

$$u(x,t)=\int_{0}^{t-\epsilon}\int_{-\infty}^{\infty} \frac{1}{\sqrt{4\pi k(t-s)}} \exp\left(-\frac{(x-y)^2}{4k(t-s)}\right)f(y,s)\,dy\,ds$$

G[x_, t_] := Exp[-x^2/(4 k t)]/Sqrt[4 π k t]

u[x_, t_] := Integrate[G[x - y, t - s] f[y, s], {s, 0, t - ϵ}, {y, -Infinity, Infinity}]

Substitute it back to the equation:

residual = D[u[x, t], t] - k D[u[x, t], {x, 2}] - f[x, t] // Simplify // Expand

It's easy to notice the last 2 terms cancel. We can manually delete them from the output, but here I'll do it programatically to make the answer more interesting:

residual2 = 
 With[{int = Integrate}, 
  residual //. 
   HoldPattern[
     coef1_ int[expr1_, rest_] + coef2_ int[expr2_, rest_]] :> 
    int[coef1 expr1 + coef2 expr2, rest]]

Remark

1. In v9.0.1 the following is enough:

   residual2 = With[{int = Integrate}, 
    residual /. 
     HoldPattern[coef_ int[int[expr_,rest1_],rest2_]]:>int[coef expr, rest2, rest1]]
   (*
   -f[x, t] + Integrate[(k f[y, t - ϵ])/(E^((x - y)^2/(4 k ϵ)) Sqrt[k ϵ]), 
                        {y, -Infinity, Infinity}]/(2 k Sqrt[Pi])
   *)
   

2. moveconst can be used in v9.0.1 in the following manner but it's a bit slow:

   residual2 = FullSimplify[D[u[x, t], t] - k D[u[x, t], {x, 2}] - f[x, t], 
     TransformationFunctions -> {moveconst, Automatic}]
   

We know (from e.g. the wiki) one possible defintion for Dirac delta is

$$ \delta_a(x)=\frac{1}{|a|\sqrt \pi}e^{-(x/a)^2}\ \text{as}\ a\to 0 $$

So the … Exp[-(…)^2] can be replaced with a … DiracDelta[…] when $\epsilon \to 0$.

It's a pity that Limit won't help in this case, as mentioned in the Possible Issues section of document of DiracDelta. Once again, it's not a bad idea to modify the output by hand, but I'll replace by coding to make the answer more interesting:

residual3 = 
 Assuming[{x ∈ Reals, ϵ > 0, k > 0}, 
  residual2 /. Exp[coef_ a_^2] :> DiracDelta[a]/Sqrt[-coef] Sqrt[Pi]]

(* -f[x, t] + f[x, t - ϵ] *)

residual3 /. ϵ -> 0
(* 0 *)

As we can see, all terms in the equation cancel, the solution is verified. (Verification for the initial condition is trivial. )

Tested in v9.0.1, v12.0.1, v12.1.0.

Answer 3

k must be both real and positive for the given domain. f[x,t] has to be well foldable as well on the domain Reals x positive Reals.

See convolve for the corresponding built-in in Mathematica.

Some theory can be found: convolution.

This book offers some more detail: Introduction to Partial Differential Equations for Scientists and Engineers using Mathematica. It is from 2014 and so fairly modern. It offers a convolution theorem in chapter 2 Integral transforms.

Handout for convolution on heat equations from Stanford University.