# How to use NDSolve with discontinuities at internal boundaries?

I don’t know how to impose discontinuous internal boundary conditions (BCs) in NDSolve, so I’ve set up an example problem to illustrate my issue. Consider the simple first-order ODE for $f(z)$ on the interval $-1 ≤ z ≤ 1$:

k f[z] + f'[z] == 0


where $k$ is an eigenvalue to be determined. At the interval edges I fix $f(-1)$ and $f(1)$ to the two constants $f_L$ and $f_R$, respectively. But I also want to be able to specify a jump discontinuity internally at $z = 0$ by setting the ratio $f(0+)/f(0-)$ to a third fixed value $r_0$. Since the above ODE is satisfied by a single exponential, it’s easy to match all the BCs and solve analytically for $k$ and $f(z)$, resulting in:

kEV[fL_, fR_, r0_] := Log[r0 fL/fR]/2;

fSol[fL_, fR_, r0_, z_] :=
Module[{k}, k = kEV[fL, fR, r0];
UnitStep[-z] fL Exp[-k (1 + z)] + UnitStep[z] fR Exp[k (1 - z)] ]


Let’s look at two examples of this solution with $f_L = e^2$ and $f_R = 1$. First, setting $r_0 = 1$ eliminates the internal discontinuity and gives a nice smooth curve with $k = 1$:

Plot[fSol[Exp[2.], 1., 1., z], {z, -1, 1},
PlotRange -> {{-1, 1}, {0, Automatic}}, Exclusions -> None,
Epilog ->
Inset[Style["k = " <> ToString[kEV[Exp[2.], 1., 1.]],
14], {-.5, .5}]]


Second, putting $r_0 = 0.5$ gives a discontinuous solution with $k = 0.653426$:

Plot[fSol[Exp[2.], 1., .5, z], {z, -1, 1},
PlotRange -> {{-1, 1}, {0, Automatic}}, Exclusions -> None,
Epilog ->
Inset[Style["k = " <> ToString[kEV[Exp[2.], 1., .5]],
14], {-.5, .5}]]


Now the crux of my question: how to get both of these results using NDSolve? It’s straightforward to duplicate the continuous plot above by simply ignoring the internal boundary and solving numerically with just the two edge BCs:

sol = NDSolve[{k[z] f[z] + f'[z] == 0, k'[z] == 0,
f[-1] == Exp[2], f[1] == 1}, {f, k}, {z, -1, 1}];

Plot[Evaluate[f[z] /. sol], {z, -1, 1},
PlotRange -> {{-1, 1}, {0, Automatic}}, Exclusions -> None,
Epilog ->
Inset[Style["k = " <> ToString[k[1.] /. sol[[1]]], 14], {-.5, .5}]]


But at this point I’m stuck and can’t see how to repeat the discontinuous case with NDSolve. According to the documentation, the “Chasing” method can solve multipoint boundary value problems, allowing the imposition of internal BCs, but the only example shown seems to assume continuity on the internal boundaries. What I need, but lack, is a mechanism in NDSolve for specifying the solution (and its derivatives for higher-order problems) just to the left and, separately, just to the right of an internal boundary. Can anyone suggest a way to make NDSolve work for my discontinuous case? Thanks!

Working solution

One can manually implement the shooting method with ParametricNDSolveValue and FindRoot:

psol = ParametricNDSolveValue[{k f[z] + f'[z] == 0, f[-1] == Exp[2],
WhenEvent[z == 0, f[z] -> r0 f[z]]}, f, {z, -1, 1}, {k, r0}];
k0 = k /. FindRoot[psol[k, 0.5][1] == 1, {k, 1}]
(*  0.653426  *)

Plot[psol[k0, 0.5][z], {z, -1, 1}]


Original -- revealing some bugginess in NDSolve

Something like this? [Edit: As @xzczx pointed out, NDSolve does this wrong. I feel it should work, but does not. I know what it looks like it's doing (see comments), but if I figure what it actually is doing, I will update.]

Block[{r0 = 0.5},
{sol} =
NDSolve[{k[z] f[z] + f'[z] == 0, k'[z] == 0, f[-1] == Exp[2], f[1] == 1,
WhenEvent[z == 0, f[z] -> r0 f[z]]},
{f, k}, {z, -1, 1}]
];

Plot[f[z] /. sol // Evaluate, {z, -1, 1}]


For some reason, NDSolve settles on a value for k[z] very close to 1, which is the solution for the OP's first BVP without the discontinuity.

• Sadly this answer is wrong. Just compare the values at z=0 with OP's, they're apparently different. – xzczd Jun 20 '15 at 12:12
• @xzczd Thanks, I guess I was tired, because it looked right last night. It is as if NDSolve solve the BVP without the WhenEvent using the shooting method and then after solving for k, it integrates with the event. I'll look into it further. It might be a bug, because it should work. (Works if I do the shooting method manually with ParametericNDSolve.) – Michael E2 Jun 20 '15 at 13:06
• @xzczd It can also be "solved" by setting the boundary condition r0 f[1] == 1, but of course that's a klugde. Still trying to see if I can track it down definitively... – Michael E2 Jun 20 '15 at 15:02
• When I first saw that ParametricNDSolve[] was finally inplemented, my first thought was "this should make shooting easier". The equivalent code using just NDSolve[] would probably have looked messy. – J. M.'s technical difficulties Jun 21 '15 at 1:10
• One possible reason for the failure of the former approach is WhenEvent can only deal with IVP (at least now). (Just checked the document of WhenEvent and saw no example for a BVP. ) – xzczd Jun 21 '15 at 4:41

You can make a change of variable to solve the problem. Here I'll use dChange for this task:

r0 = 0.5;
eqn = {k[z] f[z] + f'[z] == 0, k'[z] == 0, f[-1] == Exp[2], f[1] == 1};
c = Piecewise[{{r0, z > 0}}, 1];
neweqn = dChange[eqn, f[z] == c g[z]];
{solg, solk} = NDSolveValue[neweqn, {g, k}, {z, -1, 1}];
Plot[c solg[z], {z, -1, 1}]


• Thanks @xzczd for showing me this clever approach using dChange. Alas, I'm not sure how to use this for my real problem, where I'm trying to solve for f(z) in a fourth-order ODE with four involved internal boundary conditions interrelating f, f', f'' and f''' in a discontinuous manner. – Dr. Know Jun 20 '15 at 16:26
• @Dr.Know To deal with discontinuities of derivatives, you can try changing the independent variable, for example $f′(0+)/f′(0−)=\lambda$ can be "eliminated" by $ζ=c z$ where $c$ is a piecewise coefficient. – xzczd Jun 21 '15 at 5:44