NDSolve example - analytical vs numerical solution. How to specify initial conditions?

Question

So, I have the following equation:

y[t]^-3 == w^2 y[t]/t^2 + y''[t]

where

w^2 = (t0 k0)^2.

Analytic solution is easy to obtain:

y[t] = a Sqrt[t]

where

a = 1/Sqrt[t0 k0]

Since the solution is singular in 0, am I not sure how to specify the initial conditions for NDSolve. When I put this in NDSolve I expect to get a solution that behaves like ~t^(1/2)), but I get a ~t solution which doesn't make sense.

NDSolve[{-y''[t] - y[t]/t^2 + y[t]^(-3) == 
                0, y[?] == ?, y'[?] == ?}, y, {t, ?, ?}]

I don't have any constrains for t0 and k0 from physics (yet), I just want to understand better how the NDSolve really works.

Edit:

When I try with: w=1, k0, t0=1, it works fine:

NDSolve[{-y''[t] - 1(y[t])/(t^2) + (y[t])^(-3) == 0, y[1] == 1, y'[1] == 1},
y, {t, 1, 100}]-> ~t^1/2

behaviour.

But when I try with: w=10, k0=1, t0=10

NDSolve[{-y''[t] - 10*(y[t])/(t^2) + (y[t])^(-3) == 0, y[10] == 10, 
y'[10] == 0.045}, y, {t, 10, 100}],

I get some kind of crazy periodic functions. Why?

(IC for y'[10] == 0.045 follows directly from derivative.)

Answer 1

Mathematical overview

You seem to be using NDSolve perfectly well, and this differential equation does not seem to pose any special problems for NDSolve. There seems to be two issues raised, the oscillatory solutions and the square root solutions.

It is not uncommon for differential equations to have an isolated non-oscillatory solution among oscillatory ones. With nonlinear equations it is even less clear what the "normal" situation might be. However, we can make some rough guess at how this ODE works. We might view the differential equation in the form

$$(1)\ \ y''(t) + \left(\frac{w^2}{t^2}-\frac{1}{y(t)^4}\right)\;y(t)=0, \quad{\rm or}\quad (2)\ \ y''(t) + \frac{w^2}{t^2}\,y(t)=\frac{1}{y(t)^3}\,,$$

as a perturbation of the equidimensional equation

$$y''(t) + \frac{w^2}{t^2}\;y(t)=0\,,$$

whose solution is $y(t)=A \sqrt{t} \sin \left({\sqrt{4w^2-1} \over 2}\;\log(t) - t_0\right)$ if $w>1/2$. This in turn can be related to both $y''+ε\,y=0$ and $y''+y/(4t^2)=0$ (i.e. $w=1/2$), The first has periodic solutions if $ε>0$ and exponential growth (generically) if $ε<0$. The second has $y=a\sqrt{t}$ as a solution. If in (1) the coefficient $w^2t^{-2}-y^{-4}$ is greater than $ε>0$ over a long enough interval, then the solution must eventually decrease; and if it is less than $ε<0$ over a long enough interval, the solution must eventually increase. Further if the coefficient can remain close to $1/(4t^2)$, the solution might remain close to a square root or a logarithm times a square root.

When $y(t)$ is large, then the term $-1/y(t)^4$ becomes negligible and the coefficient will be positive. That leads to $y(t)$ decreasing. As it does, the term $-1/y(t)^4$ becomes dominant and the coefficient becomes negative.. Eventually, $y(t)$ must begin to increase. Thus oscillatory behavior seems likely to be a characteristic of many solutions. The way oscillations might be avoided is if in the transition between $y(t)$ large and $y(t)$ small, the solution does not pass across from one region to another but the coefficient approaches zero asymptotically. (It is clear from NDSolve that this does not happen.)

The ODE

Let us set up the ODE this way:

odelincoeff = (w^2/(t^2) - (y[t])^(-4));
ode = y''[t] + odelincoeff y[t]

Power function solutions

The coefficients allow the possibility that there is a solution that is a power function or a sum of power functions. To look for a power function solution, we have to find a power that works, if any. Substitute an arbitrary power function a t^n into the differential formula. The condition for a such a solution is that the result must vanish.

solcondition = ode /. y -> (a #^n &) // Expand
(*
  -a n t^(-2 + n) + a n^2 t^(-2 + n) - t^(-3 n)/a^3 + a t^(-2 + n) w^2
*)

The powers must match up in order to get cancellation. That puts a condition on n that we can use to find n.

powercondition = Simplify[
  First @ Solve[Equal @@ DeleteDuplicates @ Cases[solcondition, t^_, Infinity], n,  Reals],
  t > 0]
(*
  {n -> 1/2}
*)

Having found a suitable power, we can solve for the undetermined coefficient a:

solcondition /. powercondition
(*
  -(1/(a^3 t^(3/2))) - a/(4 t^(3/2)) + (a w^2)/t^(3/2)
*)

coeffcondition = First@Select[
   Solve[0 == solcondition /. powercondition /. t -> 1, a],
   Positive[a /. # /. w -> 1] &]
(*
  {a -> Sqrt[2]/(-1 + 4 w^2)^(1/4)}
*)

Thus a particular solution is this:

powersol = Function[t,
  Evaluate[a t^n /. powercondition /. coeffcondition // ToRadicals]]
(*
  Function[t, (Sqrt[2] Sqrt[t])/(-1 + 4 w^2)^(1/4)]
*)

For this solution the terms of coefficient of y[t] match degrees and equal the ratio -y''[t]/y[t], which of course it has to.

odelincoeff /. y -> powersol // Simplify
(*
  1/(4 t^2)
*)

Note that the coefficient is postive but approaches zero. So while the solution is concave down, it does not, for a fixed [CurlyEpsilon]>0, remain greater than [CurlyEpsilon]. Therefore, unlike the case in the general discussion above, it is not necessary for it to become decreasing, which of course it does not do.

Nearby solutions

To examine what happens near the square root solution, let us first make a change of variables, the standard substitution t = e^s for an equidimensional equation to get a linear differential form plus a nonlinear term -E^(2 s)/u[s]^3.

odelog = (ode /. y -> (u[Log[#]] &) // Factor // Numerator) / u[s]^3 /.
   {Log[t] -> s, t -> E^s} // Expand
(*
  -(E^(2 s)/u[s]^3) + w^2 u[s] - u'[s] + u''[s]
*)

We can check that our power solution still works.

odelog /. {u -> (powersol[E^#] &)} // Simplify
(*
  0
*)

We can examine solution nearby the square root solution by perturbing the solution slightly by a factor (1 + ε v[s]). Replacing ε -> 0 eliminate quadratic and higher order terms in v[s]. (It turns out that the zero-order terms cancel and the ε in the linear terms is canceled, too, by Simplify.)

odeperturbed = Simplify[
   0 == (odelog /. {u -> ((1 + ε v[#]) powersol[E^#] &)} // Together // Numerator),
   ε != 0] /. ε -> 0
(*
  4 (-1 + 4 w^2) v[s] + 4 (v′′)[s] == 0
*)

DSolve[odeperturbed, v[s], s]
(*
  {{v[s] -> E^(s Sqrt[1 - 4 w^2]) C[1] + E^(-s Sqrt[1 - 4 w^2]) C[2]}}
*)

Here we convert the exponential function to trigonometric functions. We see that as a function of s = Log[t] the solution is periodic of period2 π / Sqrt[4 w^2-1], but as a function of t, the oscillations will lengthen ast` increases.

vsol = v -> 
  Function[s, 
   Evaluate[Simplify[
       v[s] /. First @ DSolve[odeperturbed, v[s], s] // ComplexExpand // ExpToTrig, 
       w > 1 && t > 0] /.
     {C[1] -> ε (Cos[s0] + I Sin[s0])/2, C[2] -> ε (Cos[s0] - I Sin[s0])/2} // Simplify]
    ]
(*
  v -> Function[s, ε Cos[s0 + s Sqrt[-1 + 4 w^2]]]
*)

We can check that v makes a good approximation. The constant and linear terms in ε vanish, so for small ε, it will be a good approximation over some interval.

Series[
 Simplify[
  odelog /. {u -> ((1 + v[#]) powersol[E^#] &)} /. vsol /. {s0 -> 0},
  w > 1/2 && s ∈ Reals],
 {ε, 0, 2}]
% /. s -> Log[t]

(*
  -((3 (E^(s/2) (-1+4 w^2)^(3/4) Cos[s Sqrt[-1+4 w^2]]^2) ε^2)/Sqrt[2]) + O[ε]^3
  -((3 (Sqrt[t] (-1+4 w^2)^(3/4) Cos[Sqrt[-1+4 w^2] Log[t]]^2) ε^2)/Sqrt[2]) + O[ε]^3
*)

Thus the solutions near the square root will oscillate back and forth around the square root graph. Pretty cool.

Initial values for the square root

Once one obtains the correct coefficient $a$, the initial values for starting at $t = t_1$ are simply $$ y(t_1) = a \sqrt{\vphantom{f} t_1}, \quad y'(t_1) = {a \over 2 \sqrt{\vphantom{f} t_1}}\,.$$ This works for any $t_1>0$, including t0 in the question.

Numerical solution

The function sol gives the solution to the specified initial value problem. The functions rtIvals and rtSol

Clear[sol, rtIvals, rtsol];

sol[t0_?NumericQ, y0_?NumericQ, yp0_?NumericQ, w0_: 1] := 
  First @ NDSolve[{ode == 0 /. w -> w0, y[t0] == y0, y'[t0] == yp0}, 
    y, {t, 10^-20, 100}];

rtIvals[t0_, k0_] := With[{w = t0 k0},
   Simplify[{a Sqrt[t0], a/(2 Sqrt[t0])} /. First @ Solve[
       (-y''[t] - w^2 (y[t])/(t^2) + (y[t])^(-3) /. y -> (a Sqrt[#] &)) == 0
        && a > 0 && 1 - 4 w^2 < 0 && t > 0, a, Reals],
    1 - 4 w^2 < 0 && t > 0]];
rtsol[{t0_, k0_}] := sol[t0, ##, t0 k0] & @@ rtIvals[t0, k0]

Examples

Here are several solutions with w = 1. We can see the nearby solution oscillate around the red square root graph.

Plot[
 Evaluate[{y[t] /. sol[1, Sqrt[2]/3^(1/4), (1 + #) 1/(Sqrt[2] 3^(1/4))] & /@ 
    Range[-2, 2]}],
 {t, 0, 50}, 
 PlotStyle -> {Hue[0.7, 0.6, 0.6], Hue[0.5, 0.6, 0.6], Hue[0., 0.6, 1]}]

Given the nature of our analysis, a log-log plot seems like a good idea. Another approach is to subtract the square root and scale by 1/Sqrt[t], and graph as a function of s (t = E^s). The period for w = is 2 Pi/Sqrt[4 w^2 - 1]=2Pi/Sqrt[3]=3.6276. Look how well it matches. Pretty neat.

LogLogPlot[
 Evaluate[{y[t] /. sol[1, Sqrt[2]/3^(1/4), (1 + #) 1/(Sqrt[2] 3^(1/4))] & /@ 
    Range[-2, 2]}],
 {t, 1/50, 50},
 PlotStyle -> {Hue[0.7, 0.6, 0.6], Hue[0.5, 0.6, 0.6], Hue[0., 0.6, 1]}]

Plot[Evaluate[(-(4/3)^(1/4) Sqrt[E^s] + y[E^s])/Sqrt[E^s] /. 
     sol[1, Sqrt[2]/3^(1/4), (1 + #) 1/(Sqrt[2] 3^(1/4))] & /@ 
   Range[-2, 2]], {s, -Log@50, Log@50}, 
 PlotStyle -> {Hue[0.7, 0.6, 0.6], Hue[0.5, 0.6, 0.6], 
   Hue[0., 0.6, 1]}, GridLines -> {2. Pi/Sqrt[3] {-1, 0, 1}, None}]

Here are a few square root solutions for different values of t0, with k0 = 1.

Plot[Evaluate[{y[t] /. rtsol[{#, 1}] & /@ {1, 3, 9, 27}}],
 {t, 0, 50},
  PlotStyle -> {Hue[0.7, 0.6, 0.6], Hue[0.5, 0.6, 0.6], 
   Hue[0.2, 0.6, 0.6], Hue[0., 0.6, 1]}]