Using ParametricNDSolve does not work for certain parameters

Question

I want to solve the following second-order DE $$ \ddot{w}(t) H(w(t))+\frac{1}{2}(\dot{w}(t))^2 H'(w(t))+g(w(t))(c_0 (\overline t-\underline t)-w(t)) =0, $$ with boundary conditions $\dot{w}(\overline t) = 0$ and $\dot{w}(\underline t) = \frac{\overline t-\underline t}{1-G(w(\underline t))}$, and where $H(w)= G(w)(1-G(w))$, $h(w)=H'(w)$, $g(w)=G'(w)$, and $G(\cdot)$ is a cdf.

To do so, I use ParametricNDSolve as follows:

c0 = 1;
tinf = 4;
tsup = 5;

(*Uniform distribution fo G(.)*)
csup = 3; 
G[c_] = c/csup;
g[c_] = D[G[c], c];

H[w_] = G[w] (1 - G[w]);
h[w_] = D[H[w], w];

solwa = ParametricNDSolve[{w''[t]*H[w[t]] + 1/2 (w'[t])^2 h[w[t]] + g[w[t]] (c0*(tsup-tinf) - w[t]) == 0, w'[tsup] == 0, w'[tinf] == a}, {w}, {t, tinf, tsup},a]

Then, I obtain parameter $a$ as follows:

Wat = w /. solwa;
Watinf = w[a][tinf] /. solwa;
aopt = a /. FindRoot[a-(tsup - tinf)/(1 - G[Watinf]), {a,0.5}]

With the parameters given above, this works fine ($a \approx 1.1214$) and I get a solution which looks like as expected.

My problem is that, as soon as I change too much the parameters, there seems to have no solutions. For example, if I take csup=4 instead of csup=3, then I get the following error message (at the second stage when solving for $a$)

Power::infy: Infinite expression 1/0. encountered.

Power::infy: Infinite expression 1/0.^2 encountered.

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.

Power::infy: Infinite expression 1/0. encountered.

General::stop: Further output of Power::infy will be suppressed during this calculation.

Infinity::indet: Indeterminate expression ComplexInfinity+ComplexInfinity+ComplexInfinity encountered.

Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered.

General::stop: Further output of Infinity::indet will be suppressed during this calculation.

ParametricNDSolve::ndnum: Encountered non-numerical value for a derivative at t$214883 == 4.`.

ParametricNDSolve::ndnum: Encountered non-numerical value for a derivative at t$214883 == 4.`.

FindRoot::nlnum: The function value {1.05555 -1/(1-1/4 ParametricFunction[1,Internal`Bag[<1>],1,1,False,{<<7>>},{<<2>>}][1.05555][4])} is not a list of numbers with dimensions {1} at {a} = {1.05555}.

This is problematic since I want to study the solution for a broad set of parameters.

I have the same problem when using a different cdf (for instance $G(w)=1-\exp^{-\lambda w}$).

If someone can help me or point towards an alternative solution, that would be greatly appreciated.

Many thanks in advance

Answer 1

Here is a more compact approach that obtains w and a simultaneously.

NDSolveValue[{w''[t]*H[w[t]] + 1/2 (w'[t])^2 h[w[t]] + 
    g[w[t]] (c0*(tsup - tinf) - w[t]) == 0, w'[tsup] == 0, w'[tinf] == a[tinf], 
    a'[t] == 0, a[tinf] == (tsup - tinf)/(1 - G[w[tinf]])}, 
    {w, a[tinf]}, {t, tinf, tsup}, 
    Method -> {"Shooting", "StartingInitialConditions" -> {w[tinf] == .3}}]

which yields 1.12214 for a when csup = 3, and 1.08951 for a when csup = 4. See https://reference.wolfram.com/language/tutorial/NDSolveBVP.html#3518691 for details. Note that it may be necessary to adjust "StartingInitialConditions" -> {w[tinf] == .3} when other parameters in the calculation are changed.