NDSolve with parametric endpoint

Question

I want to solve numerically a system of ODEs of the form:

$$\ddot{y}+3\frac{\dot{b}}{b}\dot{y}=2Cy(y^2-1)$$ and $$\dot{b}^{2}=1-b^{2}+\frac{D}{2\sqrt{C}}b^{2}\left(\dot{y}^{2}-C(y^{2}-1)^{2}\right)$$ for the independent variables $y(t)$ and $b(t)$ and with the free parameters $C$ and $D$. The boundary conditions are: $$y(0)=0~,~~\dot{b}(0)=0~,~~b(\zeta)=0~,~~\dot{y}(\zeta)=0$$ where $\zeta$ is an unknown variable that should depend on $C$ and $D$. If $\zeta$ was well specified (e.g. $\zeta=1$) then I suspect that the system can be solved with "ParametricNDSolve" command for the parameters $C$ and $D$. However, I can't find a way to include the boundary condition with $\zeta$ unspecified. In particular, I would like to plot something like $\zeta$ vs $C$ for $D=const$.

Thank you.

Answer 1

I have edited my answer below extensively:

1. Simplified expansion about the point where b vanishes.
2. Improved accuracy of solution
3. Computed eta as a function of d for c = 1.

Integrating toward b[eta] = y'[eta] = 0 is challenging, because b[eta] occurs in the denominator of the first ODE. It is better to integrate toward b'[0] = y[0] = 0. Moreover, because the ODEs are autonomous, the origin of t can be shifted so that eta becomes 0 and 0 becomes -eta, which is convenient. To avoid using reserved symbols and also for convenience, D/(2 Sqrt[C]) is replaced by d and C by c.

ode = {b'[t]^2 == 1 - b[t]^2 + d b[t]^2 (y'[t]^2 - c (y[t]^2 - 1)^2), 
       y''[t] + 3 b'[t] y'[t]/b[t] == 2 c y[t] (y[t]^2 - 1)};

To integrate starting at b[0] = y'[0] = 0, where ode[[1]] is singular, the ODEs can be expanded about this point to obtain an approximate solution there:

CoefficientList[Normal@Series[#, {t, 0, 3}] /. {y'[0] -> 0, b[0] -> 0}, t] &
    /@ {Subtract @@ First[ode], implify[b[t] Subtract @@ Last[ode]]};
DeleteCases[Flatten[%], 0];
Solve[Thread[% == 0], Variables[%][[4 ;;]]];
{baprx, yaprx} = (Sum[Derivative[n][#][0] t^n/n!, {n, 0, 4}] /. 
    First[%] /. {y[0] -> y0, y'[0] -> 0, b[0] -> 0}) & /@ {b, y}

(* {-t + 1/6 t^3 (1 + c d - 2 c d y0^2 + c d y0^4), 
    y0 + 1/4 c t^2 y0 (-1 + y0^2) + 1/48 c t^4 y0 (-1 + y0^2) 
        (1 - c + c d + 3 c y0^2 - 2 c d y0^2 + c d y0^4)} *)

At the other boundary condition, where b'[t] = 0, ode[[2]] exhibits a branch point. Although the branch point is not an insurmountable problem, it is convenient to eliminate it by differentiating ode[[2]].

Flatten@Solve[Last[ode], y''[t]] /. Equal -> Rule;
Simplify[(D[Subtract @@ First[ode], t] /. %)/(2 b'[t])] == 0;
ode1 = {%, Last[ode]}
(* {b[t] (1 + c d - 2 c d y[t]^2 + c d y[t]^4 + 2 d y'[t]^2) + b''[t] == 0, 
    y''[t] + 3 b'[t] y'[t]/b[t] == 2 c y[t] (y[t]^2 - 1)}; *)

The additional boundary condition, for b'[0], introduced by this differentiation can be derived consistently from baprx. Now, the following short code solves the ODEs for given y0.

eps = -10^-8;
ic = {b[t] == baprx, y[t] == yaprx, b'[t] == D[baprx, t], 
    y'[t] == D[yaprx, t]} /. t -> eps;
solp = ParametricNDSolveValue[{ode1, ic, 
    WhenEvent[b'[t] == 0, "StopIntegration"]}, {b'[t], y[t], b[t], 
    y'[t]}, {t, -5, eps}, {c, d, y0}];

From this, the boundary conditions at b'[t] = 0 can be applied by (in the example c = 1, d = -0.48),

FindRoot[sf = solp[1, -.480, yf]; 
    Norm[sf[[1 ;; 2]] /. t -> (sf[[1]] /. t -> "Domain")[[1, 1]]], 
    {yf, .74, .78}, Evaluated -> False, MaxIterations -> 500]
s = solp[1, -.480, yf /. %];
tm = (s[[1]] /. t -> "Domain")[[1, 1]]
s /. t -> tm
Plot[s, {t, tm, eps}, LabelStyle -> {12, Bold, Black}, 
    PlotLegends -> Placed[{b'[t], y[t], b[t], y'[t]}, {.1, .22}]]

(* {yf -> 0.59085} *)
(* -2.28386 *)
(* {1.9082*10^-17, -1.02002*10^-15, 1.25426, 0.490884} *)

From this solution, nearby solutions (in d) can be obtained using nearby starting guesses for y0. It turns out, however, that solutions exist only for a small range of d values, at least for c = 1. Some experimentation was needed to obtain the limits on d, which are approximately d = -.47459185 with y0 = 0.764158.

and d = -.49992 with y0 = 0.0334762`.

Plots of eta, the negative of tm, and of y0 as a function of d summarize the results of several such calculations.

Second Set of Solutions

A second set of solutions can be obtained by replacing First[%] by Last[%] in the expression for {baprx, yaprx}; i.e.,

{baprx, yaprx} = (Sum[Derivative[n][#][0] t^n/n!, {n, 0, 4}] /. 
    Last[%] /. {y[0] -> y0, y'[0] -> 0, b[0] -> 0}) & /@ {b, y}

(* {t + 1/6 t^3 (-1 - c d + 2 c d y0^2 - c d y0^4), 
    y0 + 1/4 c t^2 y0 (-1 + y0^2) + 1/48 c t^4 y0 (-1 + y0^2) 
        (1 - c + c d + 3 c y0^2 - 2 c d y0^2 + c d y0^4)} *)

With this change, the expressions for solp, etc. remains unchanged but produce different solutions, as shown below. However, for convenience, define

fr[c_, d_, yg_] := Module[{yr}, 
    yr = FindRoot[sf = solp[c, d, yf]; Norm[sf[[1 ;; 2]] /. 
        t -> (sf[[1]] /. t -> "Domain")[[1, 1]]], 
        {yf, yg - .02, yg + .02}, Evaluated -> False];
    s = solp[c, d, yf /. yr];
    tm = (s[[1]] /. t -> "Domain")[[1, 1]];
{d, yf /. yr, tm, Norm[s[[;; 2]] /. t -> tm]}]

A sample solution for c = 3 and positive d is

fr[3, .13, .3]
Plot[s, {t, tm, eps}, LabelStyle -> {12, Bold, Black}, 
    PlotLegends -> Placed[{b'[t], y[t], b[t], y'[t]}, {.1, .22}]]

(* {0.12, 0.404346, -1.34058, 4.81741*10^-16} *)

The limits on d in this case are approximately d = -0.122 with y0 = 0.964856 and d = 0.166 with y0 = 0.0464389. Corresponding plots are

Finally, lots of eta, the negative of tm, and of y0 as a function of d summarize the results of several such calculations.