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$):

enter image description here

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]
  • 1
    $\begingroup$ 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. $\endgroup$
    – Nasser
    Jun 6, 2020 at 14:50
  • $\begingroup$ I checked that the formula is correct. $\endgroup$
    – Leponzo
    Jun 6, 2020 at 20:02
  • 1
    $\begingroup$ 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. ) $\endgroup$
    – xzczd
    Jun 7, 2020 at 1:58
  • $\begingroup$ @xzczd TIL the order of Integrate is the opposite. The documentation doesn't have a clear example for multiple integration. $\endgroup$
    – Leponzo
    Jun 7, 2020 at 3:34
  • $\begingroup$ I guess this design is for keeping consistent with Table. $\endgroup$
    – xzczd
    Jun 7, 2020 at 4:08

3 Answers 3


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

enter image description here

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 //. 
     coef1_ int[expr1_, rest_] + coef2_ int[expr2_, rest_]] :> 
    int[coef1 expr1 + coef2 expr2, rest]]

enter image description here


  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

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

enter image description here


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.


moveconst is from how to simplify symbolic integration thanks to celtschk

   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}] + 
    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}]


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

  • $\begingroup$ Thanks, but why is Mathematica unable to verify it by substitution though? $\endgroup$
    – Leponzo
    Jun 6, 2020 at 20:03
  • 4
    $\begingroup$ @Leponzo Because Leibniz integral rule doesn't apply in this case. $\endgroup$
    – xzczd
    Jun 7, 2020 at 1:29
  • 1
    $\begingroup$ @xzczd That was the answer I was looking for! $\endgroup$
    – Leponzo
    Jun 7, 2020 at 12:39
  • $\begingroup$ @Leponzo added verification $\endgroup$
    – Nasser
    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.

  • 1
    $\begingroup$ 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. $\endgroup$
    – MarcoB
    Jun 6, 2020 at 22:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.