Ggthjy

I need to solve some differential equations with a Dirac delta function. e.g. the source terms are like, $delta(x)$ and $partial_x delta(x)$.

Could I just use the Gaussian type function? e.g.

$$ delta(x-a) = frac{a_1}{sqrt{pi}}exp(-frac{(x-a)^2}{a_1}),$$ and,
$$ partial_x delta(x-a) = - frac{2}{sqrt{pi}}(x-a)exp(-frac{(x-a)^2}{a_1}) $$

where by setting e.g. $a_1=0.01$.

If the answer is yes, then how to choose the reasonable values of $a_1$ and $a_2$?

Thanks!

You can use a number of approximations. For example, the following three are quite often used.

diracDeltaLorentz[x_, eps_] := (eps/Pi)*(1/(x^2 + eps^2))

diracDeltaGauss[x_, eps_] := (1/(eps*Sqrt[Pi]))*Exp[-x^2/eps^2]

diracDeltaSin[x_, eps_] := (1/(Pi))*Sin[x/eps]/x



Plot[Evaluate[

  Table[diracDeltaLorentz[x, i], {i, {1/2, 1, 2}}]], 

 {x, -4, 4}, 

 PlotRange -> All, PlotStyle -> {Red, Green, Blue}]

Plot[Evaluate[

  Table[diracDeltaGauss[x, i], {i, {1/2, 1, 2}}]], 

  {x, -4, 4}, PlotRange -> All, 

  PlotStyle -> {Red, Green, Blue}]

Plot[Evaluate[

    Table[diracDeltaSin[x, i], {i, {1/2, 1, 2}}]], 

   {x, -4, 4}, PlotRange -> All, 

   PlotStyle -> {Red, Green, Blue}]

Of course, not all definitions are equivalent or well suited for integration with the given function. For example, let’s take

d1G = diracDeltaGauss[(x - 

     a*t^2/2) + (1/(2 a))*(v - a*t)^2, ϵ1]

and

d2G = diracDeltaGauss[(v - a*t) + a*(t - t0), ϵ2]

Then you can compute the integrals

 xIntDiracDeltaG = 

     Integrate[x*d1G*d2G, {x, -Infinity, Infinity}, 

      Assumptions -> {Element[{a, v, t}, Reals], 

                      ϵ1 > 0, ϵ2 > 0}]



(*Out: 

(E^(-((-a t0 + v)^2/ϵ2^2)) (2 a t - 

   v) v)/(2 a Sqrt[π] ϵ2)

*)

and then

xvIntDiracDeltaG = 

 Integrate[xIntDiracDeltaG, {v, -Infinity, Infinity}, 

  Assumptions -> {Element[{a, v, t, t0}, Reals], 

                  ϵ1 > 0, ϵ2 > 0}]



 (*Out: a t t0 - (2 a^2 t0^2 + ϵ2^2)/(4 a) *)

However if you define similar functions using the definitions of diracDeltaLorentz and diracDeltaSin above and then try to repeat computation of the same integrals, you will notice that the integrals do not converge. So, the particular choice of a diracDelta "model" depends on the functions you want to integrate.

If you need infinite differentiability, then I think a Gaussian is a good choice. For choosing parameters the parameterization given on https://en.wikipedia.org/wiki/Gaussian_function in terms of the width $sigma$ is easier:

g[{μ_, σ_}, x_] = E^(-((x-μ)/σ)^2/2)/(σ*Sqrt[2π]);

dg[{μ_, σ_}, x_] = D[g[{μ, σ}, x], x];

Choose $sigma$ smaller than any feature of your other functions.

Alternatively, if you only need first derivatives you could use a triangle function, this may be even simpler:

t[{μ_, σ_}, x_] = Piecewise[{{1/σ (1 + (x - μ)/σ), μ - σ < x < μ},

                             {1/σ (1 - (x - μ)/σ), μ < x < μ + σ}}];

dt[{μ_, σ_}, x_] = D[t[{μ, σ}, x], x];

Again, choose the width $sigma$ smaller than any feature of your other functions.

Usually the Dirac $delta$-function is defined with $mu=0$ in the above expressions (it is centered at the origin).

DSolve can often handle DiracDeltas directly, without resorting to numerical or functional approximations. For example:

DSolve[y'[x] + 7 y[x] == DiracDelta[x] + D[DiracDelta[x], x], y[x], x]

provides a nice closed form solution.