# Why the difference in fundamental solutions?

Finding a fundamental solution of a linear differential operator $$f''-c^2f$$, I try in 12.3.1 on Windows 10

DSolve[f''[x] - c^2*f[x] == DiracDelta[x - y], f[x], x]


{{f[x] -> E^(c x) C[1] + E^(-c x) C[2] - ( E^(-c x - c y) (-E^(2 c x) + E^(2 c y)) HeavisideTheta[x - y])/( 2 c)}}

To be sure, I also execute (following the documentation to DiracDelta, namely the "Properties" section, and Encyclopedia of Mathematics)

DSolve[f''[x] - c^2*f[x] == DiracDelta[-x + y], f[x], x]


{{f[x] -> E^(c x) C[1] + E^(-c x) C[2] + E^(-c x) (E^(2 c x) Inactive[Integrate][(E^(-c K[1]) DiracDelta[y - K[1]])/( 2 c), {K[1], 1, x}] + Inactive[Integrate][-((E^(c K[2]) DiracDelta[y - K[2]])/( 2 c)), {K[2], 1, x}])}}

Which of these different results is correct?

• The first one is the correct one , I think. The second solution is a general solution, which Mathematica couldn't evaluate further. May be because Mathematica doesn't recognize the symmetrie of DiracDelta: DiracDela[x]== DiracDelta[-x] ->True  but DiracDela[x-y]== DiracDelta[y-x] ->unevaluated  Nov 19, 2021 at 13:34
• My unsuccessful attempt is DSolve[f''[x] - c^2*f[x] == eps/Pi/(eps^2 + (x - y)^2), f[x], x] and then eps tends to zero from above. The result is E^(c x) C[1] + E^(-c x) C[2]. Nov 19, 2021 at 13:46
• The weak limit differs from the usual notion of the limit. Nov 19, 2021 at 13:54
• Your "unsuccessful attempt" gives a solution, similar to your first solution in your question. But last step eps->0 is not alowed! Nov 19, 2021 at 14:25
• Your last "personal" hint is misplaced, completely unnecessary and completely unfounded! Like many other users too in the stackexchange community I'll stop my contributions to your questions completely! Nov 20, 2021 at 9:11

## 2 Answers

To long for a comment:

Assuming c==1 the first solution follows to

sol1=DSolve[{f''[x] - f[x] == DiracDelta[x - y] }, f ,x][[1, 1]][[2, 2]] /. {C[1] -> 0, C[2] -> 0}



Your "unsuccessful attempt" gives

eps =.
F = Function[{x, y, eps}, Evaluate[DSolve[{f''[x] -  f[x] == eps/Pi/(eps^2 + (x - y)^2), f[0] == 0 }, f , x][[1, 1]][[2, 2]] /. {C[1] -> 0,  C[2] -> 0}]]


Plot of both solutions

Plot3D[Evaluate[{sol1, Boole[x > y] F[x,y, .1]}], {x, 0, 3}, {y, 0,3}, PlotRange -> All]


shows quite good approximation.

• You wrote "shows quite good approximation". Can you elaborate it? What exactly is well approximated? Nov 19, 2021 at 15:46
• Do you put eps=0.0;? The result of eps=.; is eps. Nov 19, 2021 at 15:51
• I find your "too long for acomment" very similar to that answer . Nov 19, 2021 at 15:57
• I added the plot command with eps=.1. Nov 19, 2021 at 15:57
• Thank you. Still don't understand all that [[1, 1]][[2, 2]] and see "a good approximation". Nov 19, 2021 at 15:59

I'm not sure why your second solution is missing the HeavisideTheta factor in the second solution, but it is fairly straightforward to show your first solution satisfies the differential equation.

eq=f''[x]-c^2*f[x]==DiracDelta[x-y]

DSolve[eq,f[x],x]//Flatten

(*   {f[x] -> C[1]*E^(c*x) + C[2]/E^(c*x) - (E^(-(c*x) - c*y)*(E^(2*c*y) - E^(2*c*x))*
HeavisideTheta[x - y])/(2*c)}   *)

f[x_]=f[x]/.%


Plug the result into the equation to test our solution.

eq = eq//Simplify

(*   2*E^(c*(x - y))*DiracDelta[x - y] + (E^(c*(y - x)) - E^(c*(x - y)))*
DiracDelta[x - y] + ((E^(2*c*x) - E^(2*c*y))*
Derivative[1][DiracDelta][x - y])/(E^(c*(x + y))*(2*c)) ==
DiracDelta[x - y]   *)


It looks like we do not have a match, but realizing that DiracDelta's only mean something inside an integral, we multiply each side by g[x] and integrate:

Integrate[eq[[1]]*g[x], {x, -Infinity, Infinity}] == Integrate[eq[[2]]*g[x], {x, -Infinity, Infinity}]
(*   ConditionalExpression[True, Element[y, Reals]]  *)


So the result satisfies the differential equation at least in the case of real y.

Doing the integrals by hand, I agree with Mathematica for the above. In the case of your second answer, doing the integrals by hand and assuming y lies inside the limits of integration I get the same result as the first solution except that it is missing the HeavisideTheta factor. The two solutions should be the same.

• Thank you for your work which confirms that the implementation of the $\delta$-distribution in Mathematica leaves much to be desired. Nov 30, 2021 at 8:26
• Replacing eq=f''[x]-c^2*f[x]==DiracDelta[x-y] by eq=f''[x]-c^2*f[x]==DiracDelta[y-x] and executing the rest of your code without any changes, I obtain ConditionalExpression[True, y \[Element] Reals] . Nov 30, 2021 at 10:43
• Moreover, ClearAll[f,x,y]; eq1 = f''[x] - c^2*f[x] == DiracDelta[x - y]; eq2 = f''[x] - c^2*f[x] == DiracDelta[y - x]; DSolve[eq2, f[x], x] // Flatten; f[x_] = f[x] /. %;eq1 = eq1 // Simplify;Integrate[eq1[[1]]*g[x], {x, -Infinity, Infinity}] == Integrate[eq1[[2]]*g[x], {x, -Infinity, Infinity}] produces ConditionalExpression[True, y \[Element] Reals], This means that the weak solution of eq2 is the weak solution of eq1. I repeat that Integrate[(E^(-c K[1]) DiracDelta[y - K[1]])/( 2 c), {K[1], 1, x}]  makes no sense in traditional math. Nov 30, 2021 at 13:19
• I really do not understand why you can't do that integral. It is in the standard form given in all the books covering DiracDelta integrals. Integrate[f[x] DiracDelta[a-x] dx. The value is f[a] if a is contained in the integration limits, and zero otherwise. For your integral the value is Exp[-c y]/(2 c) for y between 1 and x. If y is outside those values, the value is zero, and the user must actually use knowledge of the situation if y is either 1 or x. Nov 30, 2021 at 21:33