Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Symbolic solution(s) to Heat equation?

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

I am interested in variations of the heat equation:

Mathematica graphics

or more generally

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

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.

share|improve this question
4  
For those so inclined, before requesting closing this question, please provide a comment on the actual issue. Thanks. –  chris Jul 22 at 10:54
1  
While I didn't vote to close it, I suppose the reason is that you're asking why WRI does not provide functionality X. It's not something people are likely to be able to answer here. Having said that, maybe the reason they don't provide Green's function solutions is that those are convolutions and the kernel depends on the geometry, for example (eg consider the solution of the 1d heat equation in $[-1,1]$). I don't really know though. –  acl Jul 22 at 11:15
5  
@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 at 21:15
3  
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 at 21:25
2  
@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 at 8:58

3 Answers 3

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.

share|improve this answer
    
Thanks Jens that's indeed a possible generic starting point. –  chris Jul 22 at 17:34
    
Do you think we should aim at arbitrary dimensions? –  chris Jul 22 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 at 17:49
    
I am not sure I understand how the second step (to get the Green function) is general? –  chris Jul 22 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 at 20:55

Here is extensions to @Jens answer (I think) also relying on possible separation of variable. I 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]

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

operator[sol1] //FullSimplify

(* ==> 0 *)

Following exactly the same steps,

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

yields for instance:

Mathematica graphics

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]

Mathematica graphics

Heat[Δ x, 0]

Mathematica graphics

Heat[x^n, 0]

Mathematica graphics

sol1 = Heat[Δ, a, b]

Mathematica graphics

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

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

Mathematica graphics

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

Mathematica graphics

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]

Mathematica graphics

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]

Mathematica graphics

but the quadrature it returns obeys the PDE:

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

(* 0 *)

share|improve this answer
    
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 at 19:31
    
exactly : just too long for a comment. –  chris Jul 22 at 19:34
    
+1. Wow, you're doing a nice job with this. –  RunnyKine Jul 23 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 at 21:17
1  
@chris It would be easier for us readers if you would update, let's say, twice a day :) –  eldo Jul 23 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}]

Mathematica graphics

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

Mathematica graphics

We can check that this general solution satisfies the PDE

int = operator[sol1] // FullSimplify

Mathematica graphics

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

Mathematica graphics

int = operator[sol1] // FullSimplify

Mathematica graphics

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

Mathematica graphics

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

Mathematica graphics

 int = operator[sol1] // FullSimplify

Mathematica graphics

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

Mathematica graphics

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 *)

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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