Symbolic solution(s) to Heat equation?

or more generally,(eventually) Green functions to known PDEs

I am interested in variations of the heat equation:

or more generally

or even more generally (r a vector and Δ(r) a tensor)

I understand that Mathematica cannot provide a symbolic solution:

DSolve[{eqn =
D[p[x, t], t] == Δ D[p[x, t], x, x],
p[x, 0] == g[x]}, p, {x, t}]


What I do not understand is why not?

Indeed I can define a possible solution (wikipedia) as

sol = p -> Function[{x, t}, Exp[-x^2/4/t/Δ]/Sqrt[4 Pi Δ t]];


and check that it satisfies the PDE

eqn /. sol // FullSimplify

(* ==> True *)


I can even build a larger class of solution which satisfy the boundary at t=0

sol2 = p -> Function[{x, t},
Integrate[ g[y] Exp[-(x - y)^2/4/t/Δ]/Sqrt[4 Pi Δ t], {y, -Infinity, Infinity}]];


Which seems to satisfy the PDE as well, (though I don't really understand why it fails to conclude it does!)

eqn2 = eqn /. sol2 // FullSimplify


Indeed we can check by taking the limit at t-> 0 which would replace the Gaussian by a Dirac:

p[x, t] /. sol2 /. Exp[-(x - y)^2/4/t/Δ] -> DiracDelta[x - y] Sqrt[4 Pi Δ t]

(* ==> g(x) *)


Question

Could you please explain to me why no attempt is being made (by WRI or by us) along these lines?

I understand that my particular solution corresponds to a specific boundary condition, and that it might be difficult to cover all cases, but it remains surprising that this class of PDE is ignored by mathematica (e.g. those known solutions)? May be as a community we could build up a package which addresses this issue?

I truly would like to know if (i) is there indeed no general solution (ii) there is something fundamentally wrong in collecting tools providing useful (if not fully general) classes of solutions. (iii) is there something I miss which would prevent success?

Eventually, it would be great to have a Mathematica function which say would act as a lookup table and work as follow: GreenFunction[PDE,BCs] would return the corresponding Green function if it is known in the literature.

UPDATE

I am told (see comment below) that Mathematica 10.3 can now deal with the heat equation.

• For those so inclined, before requesting closing this question, please provide a comment on the actual issue. Thanks. – chris Jul 22 '14 at 10:54
• @chris One huge criticism I have with the Mathematica community is that no one wants to pick up the slack and build a package for missing functionality... everybody wants everything to come from WRI. This is not sustainable in the long run and only makes us over dependent on them. Taking this example, you think symbolic solutions to the heat equation are very important. I don't. Maybe I want a new WYSIWYG HTML and CSS editor within the Front End or maybe a useless FashionData function. We can't be blaming WRI for not aligning their business goals with an individual's personal needs can we? – rm -rf Jul 22 '14 at 21:15
• I think they brought it on themselves in part, in as much as there used to be a set of packages e.g. on mathsource which over the years they basically swallowed or sold as add ons. It certainly killed such efforts effectively. I think the fact that they include stuff and make it homogeneous with the rest of the language is great though. This site complements these past developments but does not replace them. – chris Jul 22 '14 at 21:25
• @rm-rf To make an analogy, WRI pride themselves (rightly in my opinion) to Integrate symbolically a wide range of functions. For DSolve, as far as PDE are concerned, it seems WRI writes something analogous to "we don't do integrals involving sines". You could argue that Integrals involving sines, the end users can solve for themselves, or it should be the job for the community ;-) I am trying to convince them to their DSolve function should include what is known as it would make mathematica even better. – chris Jul 23 '14 at 8:58
• DSolve can solve heat equation now in v10.3!! – xzczd Oct 13 '15 at 7:39

A first step would be to implement a convenience function that can automatically apply the method of separation of variables to separable types of equations. To show that the steps could in principle be automated, let me repeat basically the same calculation that I did for cylindrical coordinates with only slight modifications to the heat equation:

ClearAll[pt, px, x, t, p];
operator = Function[p, D[p, t] - Δ D[p, x, x]];

ansatz = pt[t] px[x];

pde2 = Expand[Apply[Subtract, operator[ansatz]/ansatz == 0]];

ptSolution =
First@DSolve[Select[pde2, D[#, x] == 0 &] == κ^2, pt[t], t];

pxSolution =
First@DSolve[Select[pde2, D[#, x] =!= 0 &] == -κ^2, px[x], x,
GeneratedParameters -> B];

ansatz /. Join[ptSolution, pxSolution]


$$C(1)\, e^{\kappa ^2 t} \left(B(1)\, e^{\frac{\kappa x}{\sqrt{\Delta }}}+B(2)\, e^{-\frac{\kappa x}{\sqrt{\Delta }}}\right)$$

The differential equation is introduced in the form operator[f] == 0, and then f is replaced by a product ansatz. The integration constants have to be named differently for the two ordinary differential equations. The separation constant is called κ. To generalize to more than two independent variables, one would also have to automate the successive introduction of integration constants, and be more careful in the identification of the terms that depend on the different variables.

Edit: Green's function

To obtain Green's function with the above starting point, one would then use the spectral representation. The eigenvalue κ is introduced blindly above, leading to an exponentially increasing time dependence. The decay factor is therefore really obtained by replacing κ by an imaginary number. But the choice in my above solution is actually more convenient in order to perform the spectral integral, because it allows me to use a trick in which the Gaussian (NormalDistribution) appears:

solution = %;

s1 =
I/σ Expectation[
solution /. {C[1] -> 1, B[1] -> 1,
B[2] -> 0}, κ \[Distributed]
NormalDistribution[k, 1/σ]] /. k -> I κ


$$\frac{i \exp \left(-\frac{x^2+2 i \sqrt{\Delta } \kappa \sigma ^2 \left(x+i \sqrt{\Delta } \kappa t\right)}{4 \Delta t-2 \Delta \sigma ^2}\right)}{\sqrt{\sigma ^2-2 t}}$$

Simplify[s1 /. σ -> 0, t > 0]


$$\frac{e^{-\frac{x^2}{4 \Delta t}}}{\sqrt{2} \sqrt{t}}$$

I didn't worry about the precise normalization factors here, just included the essential ones. What I did here is pick one of the linearly independent solutions and constructed a wave packet from it, in such a way that its limit for small width σ becomes proportional to a delta function (at $t=0$). In the Gaussian, small σ corresponds to infinite width and therefore represents the desired spectral integral. I calculate the corresponding integral using Expectation and call it s1. To check that this is also a solution (as expected from the superposition principle) you can do this:

Simplify[operator[s1] == 0]

(* ==> True *)


Then set σ to zero, to obtain the answer you found on Wikipedia.

• Thanks Jens that's indeed a possible generic starting point. – chris Jul 22 '14 at 17:34
• Do you think we should aim at arbitrary dimensions? – chris Jul 22 '14 at 17:42
• @chris Upping the dimensions is one possibility, and automating the Green's function is another. Both can be automated to some extent - which is to say that there is definitely room for Wolfram to add some functionality. But yes, it should be possible for us to do some of their work for them... with the caveat that it will probably break in less standard cases (e.g., elliptical coordinates are already much harder). – Jens Jul 22 '14 at 17:49
• I am not sure I understand how the second step (to get the Green function) is general? – chris Jul 22 '14 at 20:50
• @chris Yes, that's true, it's a trick that works for the heat equation and also the Schrödinger equation. I omitted a bunch of steps that would be needed to get there from the spectral integral as one usually finds it in textbooks. The general relation for the Green's function in terms of eigenfunctions is however also something that could be automated. Only enforcing the convergence of the integral may in general be a problem. That's where specific tricks come in handy. What I showed is the shortest possible calculation I know of, for this case (i.e., the heat eq. in the question). – Jens Jul 22 '14 at 20:55

Here is extensions to @Jens answer (I think) also relying on possible separation of variable. It is not meant as an independent answer, but complements it.

First extend his answer to 2D

ClearAll[pt, px, x, t, p];
operator = Function[p,  D[p, t] - Δ D[p, x, x] - Δ D[p, y, y]];
ansatz = pt[t] px[x] py[y];
pde2 = Expand[Apply[Subtract, operator[ansatz]/ansatz == 0]];
ptSolution =  First@DSolve[Select[pde2, (D[#, x] == 0 &&
D[#, y] == 0) &] == κ1^2 + κ2^2, pt[t], t];
pxSolution = First@DSolve[Select[pde2, D[#, x] =!= 0 &] == -κ1^2, px[x],
x, GeneratedParameters -> b1];
pySolution = First@DSolve[Select[pde2, D[#, y] =!= 0 &] == -κ2^2, py[y],
y, GeneratedParameters -> b2];

sol = ansatz /. Join[ptSolution, pxSolution, pySolution]


I can then integrate over the constants

 sol1 = Integrate[(sol /. κ1 -> I κ1), {κ1, -Infinity, Infinity}];
sol2 = Integrate[(sol1 /. κ2 -> I κ2), {κ2, -Infinity, Infinity}]


And check that this solution works

 operator[sol2] // Simplify


See also this and that solution by @Jens via separation of variable

Try Anisotropic diffusion

Clear[operator];operator[p_] := D[p, t] - Δx D[p, x, x] - Δy D[p, y, y]
ansatz = pt[t] px[x] py[y]; operator[p[t, x, y]]


(* ==> ∂p/∂t - Δx ∂^2p/∂x^2 - Δy ∂^2p/∂y^2 *)

pde2 = Expand[Apply[Subtract, operator[ansatz]/ansatz == 0]];
ptSolution = First@DSolve[Select[pde2, (D[#, x] == 0 &&D[#, y] == 0) &] ==
κ1^2 + κ2^2, pt[t], t];
pxSolution =First@DSolve[Select[pde2, D[#, x] =!= 0 &] == -κ1^2, px[x],
x, GeneratedParameters -> a];
pySolution = First@DSolve[Select[pde2, D[#, y] =!= 0 &] == -κ2^2, py[y],
y, GeneratedParameters -> b];
sol = ansatz /. Join[ptSolution, pxSolution, pySolution]
sol1 = Integrate[(sol /. κ1 -> I κ1), {κ1, -Infinity, Infinity}];
sol2 = Integrate[(sol1 /. κ2 ->I κ2), {κ2, -Infinity, Infinity}];


UPDATE

We can move to a generic coordinate system;

Let's define the Laplacian

Clear[lap];
lap[p_, coord_: "Cartesian"] :=
Laplacian[p, {x, y, z}, coord] // Expand


Let us first try and solve in Cylindrical coordinates

Clear[operator];operator[p_] := D[p, t] - Δ lap[p, "Cylindrical"]
Format[a[i_]] = Subscript[a, i]; Format[b[i_]] = Subscript[b, i];


We chose an ansatz which is mute in y (=theta) (making assumptions about the boundary condition)

ansatz = pt[t] px[x] pz[z];
pde2 = Expand[Apply[Subtract, operator[ansatz]/ansatz == 0]];
ptSolution =
First@DSolve[Select[pde2, (D[#, x] == 0 && D[#, y] == 0 &&
D[#, z] == 0) &] == κ1^2 + κ3^2, pt[t], t];
pxSolution =
First@DSolve[Select[pde2, D[#, x] =!= 0 &] == -κ1^2, px[x], x,
GeneratedParameters -> a];
pzSolution =
First@DSolve[Select[pde2, D[#, z] =!= 0 &] == -κ3^2, pz[z],
z, GeneratedParameters -> b];
sol = ansatz /. Join[ptSolution, pxSolution, pzSolution]
sol1 = Integrate[(sol /. κ1 -> I κ1), {κ1, 0, Infinity}];
sol2 = Integrate[(sol1 /. κ3 -> I κ3), {κ3, -Infinity, Infinity}]


operator[sol2] /. z -> 2 /. x -> 1 /. t -> 2 /. Δ -> 1 //N // Expand // Chop


(* 0 *)

Let's now try in spherical coordinates

Clear[operator]; operator[p_] := D[p, t] - Δ lap[p, "Spherical"]


We chose an ansatz which is mute in y,z (=theta,phi)

ansatz = pt[t] px[x] ;

pde2 = Expand[Apply[Subtract, operator[ansatz]/ansatz == 0]]
ptSolution = First@DSolve[Select[pde2, (D[#, x] == 0 && D[#, y] == 0 &&
D[#, z] == 0) &] == κ1^2, pt[t], t];
pxSolution = First@DSolve[Select[pde2, D[#, x] =!= 0 &] == -κ1^2, px[x], x,
GeneratedParameters -> a];
sol1 = Integrate[(sol /. κ1 -> I κ1), {κ1, 0,  Infinity}] // Simplify


Check that this solution is ok

operator[sol1] /. x -> 1 /. t -> 2 /. Δ -> 1 // N //  Expand // Chop

(* ==> 0 *)


Note that this works also in 2D for, e.g. Polar coordinates

Clear[operator];operator[p_] := D[p, t] - Δ Laplacian[p, {x, y}, "Polar"];
ansatz = pt[t] px[x] ;
pde2 = Expand[Apply[Subtract, operator[ansatz]/ansatz == 0]];
ptSolution = First@DSolve[Select[pde2, (D[#, x] == 0 && D[#, y] == 0 &&
D[#, z] == 0) &] == κ1^2, pt[t], t];
pxSolution =First@DSolve[Select[pde2, D[#, x] =!= 0 &] == -κ1^2, px[x], x,
GeneratedParameters -> a];
sol = ansatz /. Join[ptSolution, pxSolution];
sol1 = Integrate[(sol /. κ1 -> I κ1), {κ1, 0,Infinity}] // Simplify


operator[sol1] //FullSimplify


(* ==> 0 *)

UPDATE 2

We can move to a more general class of heat equations:

Clear[operator];
operator[p_] := D[p, t] - x Δ D[p, {x, 2}]


Note the extra x in front of Δ

ansatz = pt[t] px[x] ;
pde2 = Expand[Apply[Subtract, operator[ansatz]/ansatz == 0]]


(* ==> d pt/dt/pt(t) - (Δ x d^2px/dx^2)/ px(x) *)

 ptSolution = First@DSolve[Select[pde2, (D[#, x] == 0 && D[#, y] == 0 &&
D[#, z] == 0) &] == κ^2, pt[t], t];
pxSolution = First@DSolve[Select[pde2, D[#, x] =!= 0 &] == -κ^2, px[x], x,
GeneratedParameters -> a];
sol = ansatz /. Join[ptSolution, pxSolution];
sol1 = Integrate[(sol /. κ -> I κ), {κ, 0, Infinity}]


operator[sol1] //FullSimplify


(* ==> 0 *)

Following exactly the same steps,

operator[p_] := D[p, t] - Δ D[1/x D[p, x], x]


yields for instance:

which I think, demonstrates the potential versatility of mathematica in this context.

This can be encapsulated as a prototype of what DSolve could eventually do as follows:

Clear[Heat];
Heat[factor_: Δ, b1_: - Infinity, b2_: Infinity] :=
Module[{operator, pde2, ansatz, ptSolution, pxSolution, sol, sol1,pt, px, κ},
operator[p_] := D[p, t] - D[factor  D[p, x], x];
Print[{operator[f[t, x]] == 0, b1, b2} // TableForm];
ansatz = pt[t] px[x] ;
pde2 = Expand[Apply[Subtract, operator[ansatz]/ansatz == 0]];
ptSolution = First@DSolve[Select[pde2, (D[#, x] == 0 && D[#, y] == 0 &&
D[#, z] == 0) &] == κ^2, pt[t], t];
pxSolution = First@DSolve[Select[pde2, D[#, x] =!= 0 &] == -κ^2, px[x],
x, GeneratedParameters -> a];
sol = ansatz /. Join[ptSolution, pxSolution];
sol1 = Integrate[(sol /. κ -> I κ), {κ, b1,b2},Assumptions->t>0];
operator[sol1] /. Δ -> 1 /. x -> 2 /. t -> 3 // N //
Expand // Chop // If[# != 0, Print["not ok!"]] &; sol1];


so that, e.g.

Heat[Δ, -Infinity]


Heat[Δ x, 0]


Heat[x^n, 0]


sol1 = Heat[Δ, a, b]


soln = sol1 /. a[_] -> 1 /. C[_] -> 1 /. a -> 0 /.
b -> 1 /. Δ -> 1;

Plot[soln /. t -> 0.01, {x, -2, 2}]


ContourPlot[soln, {x, -1, 1}, {t, 0, 1}]


The anisotropic case can be encapsulated as well:

Clear[AHeat];
AHeat[factorx_: Δx, factory_: Δy, b1_: -Infinity, b2_:Infinity,b3_: -Infinity, b4_:Infinity]
:=Module[{operator, pde2, ansatz, ptSolution, pxSolution,
pySolution, sol, sol1, sol2, pt, px, py},
operator[p_] := D[p, t] - D[factorx  D[p, x], x] - D[factory  D[p, y], y];
Print[{operator[f[t, x, y]] == 0, b1, b2, b3, b4} // TableForm];
ansatz = pt[t] px[x] py[y] ;
pde2 = Expand[Apply[Subtract, operator[ansatz]/ansatz == 0]];
ptSolution = First@DSolve[Select[pde2, (D[#, x] == 0 &&
D[#, y] == 0) &] == κ1^2 + κ2^2, pt[t], t];
pxSolution = First@DSolve[Select[pde2, D[#, x] =!= 0 &] == -κ1^2, px[x],
x, GeneratedParameters -> a];
pySolution = First@DSolve[Select[pde2, D[#, y] =!= 0 &] == -κ2^2, py[y],
y, GeneratedParameters -> b];
sol = ansatz /. Join[ptSolution, pxSolution, pySolution];
sol1 = Integrate[(sol /. κ1 -> I κ1), {κ1, b1,b2}];
sol2 = Integrate[(sol1 /. κ2 -> I κ2), {κ2, b3, b4}];
operator[sol1]   /. factorx -> 1 /. factory -> 2 /. x -> 2 /.
y -> 3 /. t -> 3 // N // Expand // Chop // If[# != 0, Print["not ok!"]] &;
sol2]


so that

 AHeat[x Δx, y Δy, 0, Infinity, 0, Infinity]


UPDATE 3

Note that mathematica does provide formal solutions in cases it cannot integrate. For instance, this case has no closed form solution

 sol1 = Heat[x + x^2, 0]


but the quadrature it returns obeys the PDE:

 D[D[sol1[[1]], x] (x + x^2), x] - D[sol1[[1]], t] //Simplify// FullSimplify


(* 0 *)

• Not sure either, but I up-voted yours to cancel the downvote. Anyway, as your answer shows there is potential for generalizations. But I did say that already, too... so really this isn't a separate answer to the original question. It's more a comment. – Jens Jul 22 '14 at 19:31
• exactly : just too long for a comment. – chris Jul 22 '14 at 19:34
• +1. Wow, you're doing a nice job with this. – RunnyKine Jul 23 '14 at 20:18
• Do your comments always run this long? Just kidding. In the Heat function, you should add Assumptions->t>0 to the Integrate to speed it up and avoid ConditionalExpression in the output. – Jens Jul 23 '14 at 21:17
• @chris It would be easier for us readers if you would update, let's say, twice a day :) – eldo Jul 23 '14 at 22:05

Let me start addressing the Green function part of the question.

Lets define a Heat equation and its generic solution (see above)

operator[p_] := D[p, t] - D[Δ  D[p, x], x];
sol = Heat[Δ, -Infinity]


and build a general solution via superposition as:

sol1 = Integrate[(sol /. x -> x - y) g[y], {y, -Infinity, Infinity}]


Plot[sol1 /. g -> Function[x, Exp[-x^2/2]] /. a[_] -> 1 /.
C[_] -> 1 /. Δ -> 1 /. t -> Range[4] //   Evaluate, {x, -8, 8}]


We can check that this general solution satisfies the PDE

int = operator[sol1] // FullSimplify


Modulo a little help to mathematica for simplification:

int /. Integrate -> Int /. a_ Int[b_, c__] -> Int[a b , c] /.
Int[a_, c__] + Int[b_, c__] -> Int[a + b, c] /.
Int -> Integrate // Simplify


(* ===> 0 *)

We can also add a shift in time:

sol1 = Integrate[(sol /. x -> x - y /. t -> t - τ) g[y, τ],
{y, -Infinity, Infinity}, {τ, 0, Infinity}]


int = operator[sol1] // FullSimplify


which needs a bit of help to simplify to zero (why?)

int /. Integrate -> Int //. a_ Int[b_, c__] -> Int[a b , c] //.
Int[a_, c__] + Int[b_, c__] :>  Int[Simplify[a + b], c] /.
Int[Int[a_, b__], c__] :>  Int[a, Sequence @@ Sort[{b, c}]] //.
a_ Int[b_, c__] -> Int[a b , c] //.
Int[a_, c__] + Int[b_, c__] :>  Int[Simplify[a + b], c] /.
Int -> Integrate // Simplify


(* ===> 0 *)

It works for general classes of solution with less trivial boundary condition and diffusion coefficient

Clear[sol];
operator[p_] := D[p, t] - D[Δ x D[p, x], x];
sol = Heat[x Δ, 0] /. a[_] -> 1
sol1 = Integrate[(sol /. x -> x - y) g[y], {y, 0, Infinity}]


int = Simplify[operator[sol1] /. Integrate -> Int/.
a_ Int[b_, c__] -> Int[a b , c] //. Int[a_, c__] +
Int[b_, c__] -> Int[a + b, c]] /. Int -> Integrate // FullSimplify


Or for 2D diffusion as well:

operator[p_] :=  D[p, t] - D[Δ  D[p, x], x] - D[Δ  D[p, y], y];
sol = AHeat[ Δ, Δ];
sol1 = Integrate[(sol /. x -> x - x1 /. y -> y - y1) g[x1,
y1], {x1, -Infinity, Infinity}, {y1, -Infinity, Infinity}]


 int = operator[sol1] // FullSimplify


though it requires a bit more sweat to show it nullifies the operator

int /. Integrate -> Int //. a_ Int[b_, c__] -> Int[a b , c] //.
Int[a_, c__] + Int[b_, c__] :>  Int[Simplify[a + b], c] /.
Int[Int[a_, b__], c__] :>  Int[a, Sequence @@ Sort[{b, c}]] //.
a_ Int[b_, c__] -> Int[a b , c] //.
Int[a_, c__] + Int[b_, c__] :>  Int[Simplify[a + b], c] /.
Int -> Integrate // Simplify


(* ===> 0 *)

Finally for anisotropic diffusion with less trivial initial condition:

operator[p_] := D[p, t] - D[Δx x  D[p, x], x] -D[Δy y  D[p, y], y];

sol = AHeat[ Δx x, Δy y, 0, Infinity, 0, Infinity];
sol1 = Integrate[(sol /. x -> x - x1 /. y -> y - y1) g[x1, y1], {x1,
0, Infinity}, {y1, 0, Infinity}];
int = operator[sol1];
int2 = int /. Integrate -> Int //.
a_ Int[b_, c__] -> Int[a b , c] //.
Int[a_, c__] + Int[b_, c__] :>  Int[Simplify[a + b], c] /.
Int[Int[a_, b__], c__] :>  Int[a, Sequence @@ Sort[{b, c}]] //.
a_ Int[b_, c__] -> Int[a b , c] //.
Int[a_, c__] + Int[b_, c__] :>  Int[Simplify[a + b], c]/. Int-> Integrate // Simplify


Now Its tricky to simplify the above but we can cheat by looking at it for specific values:

int2 /. a[_] -> 1 /. b[_] -> 1 /. Δx -> 1 /. Δy -> 2 /.
g[x_, y_] -> DiracDelta[x - 2] DiracDelta[y - 1]


(* ===> 0 *)