# Find height of oscillating solution from NDSolve

There are a class of differential equations whose solution oscillates around some equilibrium indefinitely after settling down. Below is an implementation of such an equation

t0 = 10^-6; tf = 50;
y0 = 2; dy0 = 0;
sol = NDSolve[{y''[t] == -(0.1 + 1/t) y'[t] - 1 + 6*y[t]^2 - 4 y[t]^3,
y[t0] == y0, y'[t0] == dy0}, {y}, {t, t0, tf}];


which produces an oscillating solution around y = 1.3.

I'd like to do two things. First, I'd like to stop the integration once a steady-state equilibrium is achieved. Second, I would like the find the midpoint of the oscillations, i.e., the height (for this, I could simply average over the late-time solutions so it's less of a problem).

In more complicated problems, the frequency might be increasing indefinitely, but the equilibrium height remains the same. I would hope to be able to accommodate this in the final solution.

Is there a way to easily achieve this using WhenEvent or something similar?

• You can use the Hilbert transform to shift the oscillation by a half-period. Take then square root of the sum of squares to obtain the envelope, mathematica.stackexchange.com/a/49081/9469 Commented Jan 6 at 13:39

## 2 Answers

I am not sure if I understand the question--are you looking for a way to collect the maxima? If so, you could do this:

{sol, extrema} = Reap[
NDSolveValue[{y''[t] == -(1/10 + 1/t) y'[t] - 1 + 6*y[t]^2 -
4 y[t]^3, y[t0] == y0, y'[t0] == dy0,
WhenEvent[y'[t] == 0, Sow[{t, y[t]}]]}, y, {t, t0, tf}]
]


And then:

ListPlot[extrema,AxesLabel->{"x","extremal values"}]


or

ListPlot[extrema[[1, 2 ;; -1 ;; 2]],{"x","maximal values"}]

Show[
Plot[sol[t], {t, t0, tf}, PlotRange -> All,
AxesLabel -> {"x", "Maximal Values"}],
ListPlot[extrema[[1, 2 ;; -1 ;; 2]],
PlotStyle -> {AbsolutePointSize[2], Red}, PlotRange -> All]
]


The possible stationary values (y'[t]==0,y''[t]==0) from your ode follow from

asy = Values@Solve[-1 + 6*y[t]^2 - 4 y[t]^3 == 0, y[t], Reals] //Flatten
asy//N (*{0.5, -0.366025, 1.36603}*)


To realize WhenEvent try (I have increased the value of tf)

t0 = 10^-6; tf = 250;
y0 = 2; dy0 = 0;
Y = NDSolveValue[{y''[t] == -(0.1 + 1/t) y'[t] - 1 + 6*y[t]^2 -
4 y[t]^3, y[t0] == y0, y'[t0] == dy0,
WhenEvent[
RealAbs[ y'[t]] < 0.001 && RealAbs[y[t] - asy[[3]]] < .001 ,
"StopIntegration"]}, y , {t, t0, tf}];
Plot[Y[t], {t, t0, Y["Domain"][[1, 2]]}, PlotRange -> All,
GridLines -> {None, grid}]


addendum

For the case you have no apriori information about the asymptotes check condition |D[y[t]^2,t]| <<1

Try

Y = NDSolveValue[{y''[t] == -(0.1 + 1/t) y'[t] - 1 + 6*y[t]^2 -
4 y[t]^3, y[t0] == y0, y'[t0] == dy0,
WhenEvent[   RealAbs[y[t] y'[t]] < .000001 , "StopIntegration"]},
y , {t, t0, tf}];
Plot[Y[t], {t, t0, Y["Domain"][[1, 2]]}, PlotRange -> All,
GridLines -> {None, grid}]
`
• Thank you for the solution. However, I am more interested in the abstract case where one can't necessarily use the differential equation to deduce where the asymptotes are. Is it possible solve this without first knowing through some other means? Commented Jan 5 at 17:07
• How would you define asymptotes in general ? Commented Jan 5 at 17:28
• The solution goes as $f(t) = c + A(t)\cos[\omega(t) + \delta]$ as the solution goes to $t\to\infty$. The asymptotes would be $c$ in this case. Commented Jan 8 at 8:01
• @MarcosMFlores Did you recognize the addendum in my answer? Commented Jan 8 at 9:51
• I don't think the condition you specified will always work. For example, consider the function $y(t) = 2 + \cos(t^2)/(t + 1)$ for $t > 0$. The derivative of $y^2(t)$ is not small, but the function still approaches a fixed value of 2. Commented Jan 8 at 10:40