# Boundary value problem with a DiracDelta

It seems that Mathematica can solve an initial value problem with a DiracDelta, but not a boundary value problem with a DiracDelta. Is there a workaround?

Consider for example a differential equation with this form:

f''[x]+x f'[x]== x^2 DiracDelta[x-1]


in the range 0,<x<10 with the boundary conditions f[0]=0, f[10]=1

• What do you mean by can't solve? I do get a solution. Also, for a one-dimensional case, there is no difference between boundary conditions and initial conditions. – Roderic Dec 9 '20 at 19:17
• @Roderic Maybe I chose the wrong example that can be solved with DSolve, but I need a numerical solution and NDSolve does cannot solve the above equation as the image shows. – mattiav27 Dec 9 '20 at 19:24
• right DSolve does give a solution in terms of Erf and HeavisideTheta. I don't know much about the numerical algorithms, sorry. – Roderic Dec 9 '20 at 19:26
• Manual shooting method? – Michael E2 Dec 9 '20 at 19:30
• @MichaelE2 could you post an answer? I am not good with numerical methods either – mattiav27 Dec 9 '20 at 19:31

@MichaelE2 gave the idea to use a shooting method, because NDSolve is only able to handle initial value problems involving DiracDelta

First solve the problem with a parametric slope f'[0]==fs0

F = ParametricNDSolveValue[{f''[x] + x f'[x] == x^2*DiracDelta[x - 1],
f[0] == 0, f'[0] == fs0}, f , {x, 0, 10}, fs0]


Now choose fs0 to fullfill the second boundary condition f[10]==1

sol = FindRoot[F[fs0][10] == 1, {fs0, 1}]

Plot[F[fs0 /. sol][x], {x, 0, 10}, PlotRange -> All]


• aargh, I was just entering my answer when yours appeared. :) – Michael E2 Dec 9 '20 at 21:41
• @MichaelE2 Sorry – Ulrich Neumann Dec 9 '20 at 21:56
• No problem. ;-) – Michael E2 Dec 9 '20 at 21:57
• @MichaelE2 I didn't know NDSolve can handle DiracDelta when dealing with IVP of ODE! Further test shows this seems to be related to the setting Method -> "DiscontinuityProcessing" -> True, which is added since v9. – xzczd Dec 10 '20 at 3:27

Replacing DiracDelta by its approximation in the weak topology helps:

s = NDSolve[{f''[x] + x f'[x] == x^2 *0.01/Pi/((x - 1)^2 + 0.01^2),
f[0] == 0, f[10] == 1}, f[x], {x, 0, 10}]
Plot[f[x] /. s, {x, 0, 10}, PlotRange -> All]


• Rigorously saying, we replace the δ-distribution which is not associated with any usual function by the distribution associated with a usual function 0.01/Pi/((x - 1)^2 + 0.01^2). – user64494 Dec 9 '20 at 20:01
ClearAll[f, x];
psol = ParametricNDSolveValue[{f''[x] + x f'[x] ==
x^2 DiracDelta[x - 1], f[0] == 0, f'[0] == p},
f, {x, 0, 10}, {p}];

FindRoot[psol[p][10] == 1, {p, 0.2}]
bvpsol = psol[p] /. %;
(*  {p -> 0.274728}  *)

bvpsol[{0, 10}]
% - {0, 1}
(*
{0., 1.}
{0., -1.44329*10^-15}
*)