# Inhomogeneous Heat Equation Formula Not Satisfying IVP

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]
• btw, you need to also add k>0 if you are going to do it this way, but this does not help with the integration part there. also on V 12.1 I get no Indeterminate in the result, just unevaluated integrals. You might want to double check the Green function you used again just in case. Commented Jun 6, 2020 at 14:50
• I checked that the formula is correct. Commented Jun 6, 2020 at 20:02
• The definition of u should be Integrate[G[x - y, t - s]*f[y, s], {s, 0, t}, {y, -Infinity, Infinity}] We cannot arbitrarily interchange the order of integration in this case. (This doesn't help solving the problem, though. ) Commented Jun 7, 2020 at 1:58
• @xzczd TIL the order of Integrate is the opposite. The documentation doesn't have a clear example for multiple integration. Commented Jun 7, 2020 at 3:34
• I guess this design is for keeping consistent with Table. Commented Jun 7, 2020 at 4:08

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.

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

• Thanks, but why is Mathematica unable to verify it by substitution though? Commented Jun 6, 2020 at 20:03
• @Leponzo Because Leibniz integral rule doesn't apply in this case. Commented Jun 7, 2020 at 1:29
• @xzczd That was the answer I was looking for! Commented Jun 7, 2020 at 12:39
• @Leponzo added verification Commented Jun 7, 2020 at 14:08

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.

• Can you show what you mean by using convolve? It would be best if you included code. Link-only answers become useless when the links invariably break at some point in the future. Commented Jun 6, 2020 at 22:04