# Syntax Problem In NDSolve with Free Boundaries

I am considering a problem with two free boundaries. When I consider the equation:

ESG2[mL0_?NumberQ, mH0_?NumberQ] :=
Module[{mL = mL0, mH = mH0}, Clear[v, f, d];
f = v[m] /. (NDSolve[{v''[m] + v'[m] + v[m] + 0.01 mH == 0,
Derivative[1][v][0] == 0.6 , v[mH] == v[mL]*1.01},
v[m], {m, 0, mH},
Method -> {"Shooting",
"StartingInitialConditions" -> {v'[0] == 1, v[0] == 8}}][[
1]]);
dL = D[f, m] /. {m -> mL};
dH = D[f, m] /. {m -> mH}; {dL - 1/2, dH } ]
sol = FindRoot[ESG2[mL, mH], {{mL, .3}, {mH, 1}}]



Then the program solves it beautifully. However, apparently this is not the right syntax when the equation to solve specifies the value function evaluated in one of the free boundaries, namely:

ESG2[mL0_?NumberQ, mH0_?NumberQ] :=
Module[{mL = mL0, mH = mH0}, Clear[v, f, d];
f = v[m] /. (NDSolve[{v''[m] + v'[m]+ v[m] + v[mL] + 0.01 mH == 0,
Derivative[1][v][0] == 0.6 , v[mH] == v[mL]*1.01},
v[m], {m, 0, mH},
Method -> {"Shooting",
"StartingInitialConditions" -> {v'[0] == 1, v[0] == 8}}][[
1]]);
dL = D[f, m] /. {m -> mL};
dH = D[f, m] /. {m -> mH}; {dL - 1/2, dH } ]
sol = FindRoot[ESG2[mL, mH], {{mL, .3}, {mH, 1}}]



In this second case, the error message reads:

NDSolveShooting::idelay: -- Message text not found --
ReplaceAll::reps: {0.01 +v[0.3]+(v^\[Prime])[m]+(v^\[Prime]\[Prime])[m]==0,(v^\[Prime])[0]==0.6,v[1.]==1.01 v[0.3]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.


Any suggestion of how to modify the syntax of the second problem to make it work is greatly appreciated!

• NDSolve seems to think that your ODE is a delay differential equation. Work around this with DSolveValue[{v''[m] + v'[m] + v[m] + c == 0, Derivative[1][v][0] == 6/10, v[mH] == v[mL]*101/100}, v[m], m] // FullSimplify and then substitute for c. Commented Mar 3, 2022 at 19:29

As noted in a comment, the second code block in the question fails because NDSolve is not capable of solving the ODE, because v[mL] is needed to begin integration at m = 0. However, since the ODE is linear, DSolve can, with a little help, obtain a symbolic solution for either code block. Begin by demonstrating this approach for the first code block.

sc = DSolveValue[{v''[m] + v'[m] + v[m] + mH/100 == 0, v'[0] == 6/10,
v[mH] == v[mL]*101/100}, {v'[mL] - 1/2, v'[mH]}, m] // FullSimplify;


the result of which is a bit long to show here. Numerical roots are quickly obtained from

FindRoot[sc, {{mL, .3}, {mH, 1}}]
(* {mL -> 1.81437, mH -> 6.73444} *)


as in the question. However, this is by no means the only root, as can be seen from

ContourPlot[{sc[[1]] == 0, sc[[2]] == 0}, {mL, 0, 8}, {mH, 0, 16},
RegionFunction -> Function[{mL, mH}, mL < mH], PlotPoints -> 100,
FrameLabel -> {mL, mH}, LabelStyle -> {15, Bold, Black}]


Every intersection is a root of sc. Now, it turns out that DSolve cannot solve the modified ODE in the second block of code, but it can solve

DSolveValue[{v''[m] + v'[m] + v[m] + c + mH/100 == 0, v'[0] == 6/10,
v[mH] == v[mL]*101/100}, {v'[mL] - 1/2, v'[mH], v[mL]}, m] // FullSimplify;


from which c == v[mL] gives the desired results.

sc = Simplify[Most[%] /. Flatten[Solve[c == Last[%], c]]];

FindRoot[sc, {{mL, .3}, {mH, 1}}]
(* {mL -> 1.81254, mH -> 6.73368} *)


which is nearly identical to the corresponding root for the first ODE. In fact, plotting sc for the second ODE yields curves indistinguishable for those above. Why this should be so follows from the following short analysis.

{v'[mL] - 1/2, v'[mH]} /. v -> Function[m, w[m] - c]
{v''[m] + v'[m] + v[m] + c + mH/100 == 0, v'[0] == 6/10,
v[mH] == v[mL]*101/100} /. v -> Function[m, w[m] - c]

(* {-(1/2) + w'[mL], w'[mH]}
{mH/100 + w[m] + w'[m] + w''[m] == 0, w'[0] == 3/5,
-c + w[mH] == 101/100 (-c + w[mL])} *)


In other words, the transformation v -> Function[m, w[m] - c] eliminates c from the second code block, except for the last expression, -c + w[mH] == 101/100 (-c + w[mL]). But, since c is small compared to 100 w[mL], c is almost transformed away.

I am confident that such ODEs can be solved even if nonlinear,but doing so is beyond the scope of this answer.

• Thanks a lot bbgodfrey for this. Unfortunately, the ODE I have is actually more complex that I put here-- it was very long and it is non linear (as you comment in your last bit), which is why for simplicity I did not put it here. So I need to solve the problem numerically-- Mathematica cannot find an analytical solution for my ODE in the first place. If you have any further insights on your last comment (which is precisely what I have to do!!) I m deeply grateful. Ultimately I should do a transformation like the one you suggest, but I cannot use DSolveValue.
– MCB
Commented Mar 4, 2022 at 14:24
• Using your logic, I solved it with a new variable x ESG2[mL0_?NumberQ, mH0_?NumberQ, x0_?NumberQ] := Module[{mL = mL0, mH = mH0, x = x0}, Clear[v, f, d]; f = v[m] /. (NDSolve[{v''[m] + v'[m] + v[m] + x + 0.01 mH - 0.01 mL == 0, Derivative[1][v][0] == 0.6 , v'[mH] == 0}, v[m], {m, 0, mH}, Method -> {"Shooting", "StartingInitialConditions" -> {v'[0] == 1, v[0] == 8}}][[ 1]]); dL = D[f, m] /. {m -> mL}; fH = f /. {m -> mH}; fL = f /. {m -> mL}; {dL - 1/2, fH - x*2, x - fL} ] sol = FindRoot[ESG2[mL, mH, x], {{mL, .3}, {mH, 1}, {x, 9}}]
– MCB
Commented Mar 4, 2022 at 14:40
• @MCB Yes, this is the approach I had in mind in the last sentence of my answer above. I suggest that you post it as an answer. I noticed that you have replaced dH by fH - x*2`, which differs from the second block of code in your question. Commented Mar 5, 2022 at 18:25