# NDSolve with parametric endpoint

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.

• Edit your question to include a specific example of the problem rather than a vague description. Mar 13 at 4:53
• Probably, solutions exist only for some C and D. For what ranges of these parameters do you expect solutions? Mar 14 at 1:51

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' = y = 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 = y' = 0, where ode[] is singular, the ODEs can be expanded about this point to obtain an approximate solution there:

CoefficientList[Normal@Series[#, {t, 0, 3}] /. {y' -> 0, b -> 0}, t] &
/@ {Subtract @@ First[ode], implify[b[t] Subtract @@ Last[ode]]};
DeleteCases[Flatten[%], 0];
{baprx, yaprx} = (Sum[Derivative[n][#] t^n/n!, {n, 0, 4}] /.
First[%] /. {y -> y0, y' -> 0, b -> 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[] exhibits a branch point. Although the branch point is not an insurmountable problem, it is convenient to eliminate it by differentiating ode[].

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', 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[] /. t -> "Domain")[[1, 1]]],
{yf, .74, .78}, Evaluated -> False, MaxIterations -> 500]
s = solp[1, -.480, yf /. %];
tm = (s[] /. 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][#] t^n/n!, {n, 0, 4}] /.
Last[%] /. {y -> y0, y' -> 0, b -> 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[] /. t -> "Domain")[[1, 1]]],
{yf, yg - .02, yg + .02}, Evaluated -> False];
s = solp[c, d, yf /. yr];
tm = (s[] /. 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.  • Thank you for this code! It seems to work pretty well but as you say it can be better optimized. If I figure something better I will post. I do think there is a missing factor of y0 in the expression for y'[eps]. Also, when I run the code the variable tm does not acquire a value and quite frankly I don't know how to extract it. Could you please clarify that? Finally, I am pretty sure that this can be expanded with ParametricNDSolve somehow. Mar 15 at 3:28
• @GeorgeFanaras There was a missing factor of y0 in the expression for y'[eps]. I have corrected this and made the necessary minor changes in the sample result. Also, it is straightforward to use ParametricNDSolve, and doing so would make good sense in a more automated version of my answer. I may try to produce an automated version, when I find time. Mar 15 at 10:18
• Actually the algorithm you wrote is incomplete since you are fixing the value y0. However, this value does not necessarily satisfy the boundary conditions. The correct approach would be write WhenEvent[Abs[Re[b'[t]]] < eps && Abs[Re[y[t]]]<eps , tm = t; "StopIntegration"] and use a shooting method that finds the value of y0 that satisfies it. I have not been able to figure out how to write this yet. Mar 16 at 2:37
• @GeorgeFanaras Actually what I did for the sample result in my answer was to choose y0, find the value of tm at which y = 0, and then adjust d to cause b' = 0 there. Not optimal, but I ran out of time to automate the process. Mar 16 at 3:20
• @GeorgeFanaras My update partially automates the solutions and obtains eta vs d for c = 1`. I hope that you find this approach helpful. Mar 28 at 1:49