# Solve PDE with DiracDelta function

This question is related to Analytic solution of dynamic Euler–Bernoulli beam equation with compatibility condition. I think it is more appropriate to open another question on this topic.

In the question, the compatibility conditions can be introduce into equations using DiracDelta. For example, for the spring compatibility condition, the governing equation can be rewritten as

$$EI\cfrac{\partial^4 w}{\partial x^4} + \mu\cfrac{\partial^2 w}{\partial t^2} + k \delta(x-L/2) w(x,t) = 0$$

Then, we try to use Laplace Transform mentioned by @xzczd,

eqn = EI D[y[x, t], {x, 4}] + μ D[y[x, t], {t, 2}] +
k DiracDelta[x - L/2] y[x, t] == 0;
ic = {y[x, 0] == Sin[x/L Pi], Derivative[0, 1][y][x, 0] == 0};
bc = {y[0, t] == 0, y[L, t] == 0, Derivative[2, 0][y][0, t] == 0,
Derivative[2, 0][y][L, t] == 0};
teqn = With[{l = LaplaceTransform},
l[{eqn, bc}, t, s] /. HoldPattern@l[u_, t, s] :> u] /. Rule @@@ ic


This step goes well, the equation is transformed to s-domain, but with DiracDelta function.

Apply DSolve to the equation returns 0

DSolve[teqn, y[x, t], x][[1, 1, -1]]


So, I thought maybe use FourierTransform to transform $x$ to $\omega$ domain may solve this problem. However, Mathematica do not return a result from FourierTransform

FourierTransform[eqn, x, s]


There are several questions related to this question but they are mostly using NDSolve. Here are two examples:

Solving a PDE containing DiracDelta

How to correctly use DSolve when the force is an impulse (dirac delta) and initial conditions are not zero

• I don't think using the delta function is practical here. It's probably better to use the general solution in each domain (left and right half) and match them according to the compatibility conditions implied by the delta function. That should reduce the problem to a set of simultaneous (algebraic, not differential) equations for the integration constants in each interval once you know the general solution as derived in my answer to the linked question.
– Jens
Commented May 26, 2015 at 4:48
• In v9.0.1, DSolve returns unevaluated after a warning. For the FourierTransform part: unlike LaplaceTransform, currently FourierTransform isn't handy. To circumvent this, here's a shell for enhancing. But what's really troublesome is: the ODE isn't defined in a infinite domain, so Fourier transform can't be used directly, theoretically one can extend the domain in cycles with e.g. Mod, but personally I never succeeded in the subsequent step: perhaps a artificial periodic function is too hard for FourierTransform, I'm not sure. Commented May 26, 2015 at 5:20
• @Jens substituting boundary condition into the general solution and solve the equations analytically in Mathematica is not a easy task to do. Actually, it will be very soon end up to an eigenvalue problem for this problem. Commented May 26, 2015 at 5:33
• @xzczd it seems this problem has to fall back to eigenvalue problem. And the best we can get is a not very long eigenvalue equation. Commented May 26, 2015 at 5:45
• Is the product of one-dimensional delta distribution by a function of two variables DiracDelta[x - L/2]*y[x, t] defined at all? Commented Jan 3, 2020 at 19:08

This is not a total answer, but it may get you a little farther with the FourierTransform. The pde:

eqn = EI D[y[x, t], {x, 4}] + μ D[y[x, t], {t, 2}] +
k DiracDelta[x - L/2] y[x, t] == 0


Although MMa doesn't like taking the ft of the entire pde at once, we can do it term by term.

fteq = EI FourierTransform[D[y[x, t], {x, 4}], x, s] + μ FourierTransform[D[y[x, t], {t, 2}], x, s] +
k FourierTransform[y[x, t] DiracDelta[x - L/2], x, s] == 0
(*EI*s^4*FourierTransform[y[x, t], x, s] +
μ*FourierTransform[Derivative[0, 2][y][x, t], x, s] +
(k*E^((I*L*s)/2)*y[L/2, t])/Sqrt[2*Pi] == 0*)


Converting this by hand to an ode MMa can work with:

fteq = EI s^4 ys[t] + μ ys''[t] + y[L/2, t] (k E^((I L s)/2))/Sqrt[2 π] == 0

sol = DSolve[fteq, ys[t], t] // Flatten


We get a solution, but not knowing y[L/2,t] this is about as far as we can go. MMa can invert the FourierTransform for a constant value for y[L/2,t] and presumably for other functions of t, but with the DiracDelta function at x = L/2, it is hard to know the value of y there.

• Hmm… I'm afraid this answer isn't correct, because FourierTransform is a tool for solving differential equation defined in a unbounded domain, while OP's problem is bounded in $[0,L]$. It might be possible to first extend the domain periodically and then solve the problem with FourierTransform, I myself never succeeded though… Commented Aug 12, 2018 at 1:36

After

eqn = EI D[y[x, t], {x, 4}] + μ D[y[x, t], {t, 2}] + k DiracDelta[x - L/2] y[x, t] == 0;
ic = {y[x, 0] == Sin[x/L Pi], Derivative[0, 1][y][x, 0] == 0};
bc = {y[0, t] == 0, y[L, t] == 0, Derivative[2, 0][y][0, t] == 0, Derivative[2, 0][y][L, t] == 0};
teqn = With[{l = LaplaceTransform}, l[{eqn, bc}, t, s] /. HoldPattern@l[u_, t, s] :> u] /. Rule @@@ ic


we need to do

teqns = teqn /. {t-> s}


and then

DSolve[teqns, y, {x, s}]


To have a result, suppose that instead we have (without the DiracDelta)

teqns = \[Mu] (-s Sin[(\[Pi] x)/L] + s^2 y[x, s]) + EI D[y[x,s],{x,4}] == 0


so following with

yxs = y[x,s] /. DSolve[teqns, y, {x, s}]


we obtain an output that should be concluded by

yxt = InverseLaplaceTransform[yxs, s, t]


There are many issues to handle. It seems that DSolve doesn't like the DiracDelta in the PDE and also the boundary conditions are not so easy to handle.

• DSolve[teqns, y, {x, s}] returns the input in version 12.0 on Windows 10. Commented Aug 31, 2020 at 5:19
• It gives to me with the same version $$\left\{c_1(s) e^{\frac{(-1)^{3/4} \sqrt[4]{\mu } \sqrt{s} x}{\sqrt[4]{\text{EI}}}}+c_2(s) e^{-\frac{\sqrt[4]{-1} \sqrt[4]{\mu } \sqrt{s} x}{\sqrt[4]{\text{EI}}}}+c_3(s) e^{-\frac{(-1)^{3/4} \sqrt[4]{\mu } \sqrt{s} x}{\sqrt[4]{\text{EI}}}}+c_4(s) e^{\frac{\sqrt[4]{-1} \sqrt[4]{\mu } \sqrt{s} x}{\sqrt[4]{\text{EI}}}}+\frac{\mu L^4 s \sin \left(\frac{\pi x}{L}\right)}{\left(\pi ^2 \sqrt{\text{EI}}-i \sqrt{\mu } L^2 s\right) \left(\pi ^2 \sqrt{\text{EI}}+i \sqrt{\mu } L^2 s\right)}\right\}$$ Commented Aug 31, 2020 at 6:53
• Now if $c_1(s)=c_2(s)=c_3(s)=c_4(s)=0$ due to boundary conditions we will get $$y(x,t)=\sin \left(\frac{\pi x}{L}\right) \cos \left(\frac{\pi ^2 \sqrt{\text{EI}} t}{\sqrt{\mu } L^2}\right)$$ Commented Aug 31, 2020 at 7:06