NDSolve doesn't work even though DSolve does

Question

I'm trying to solve a differential equation. Using the following code

Eq=A'[t]^2==1/A[t]-1
sol=DSolve[{Eq, A[0] == 1}, A[t], {t, 0, 2}]
Plot[A[t] /. sol[[1]], {t, 0, 2}]

gives the correct plot of the solution. The reason I write sol[[1]] is because DSolve outputs 2 solutions, one for t>0 and the other for t<0, and I'm only interested in the first one.

However, if I use NDSolve instead, the two solutions given are just straight lines, A[t]=1, which is technically a solution but not the solution I'm looking for and not the solution DSolve gives, which is a curve. What's happening here and what can I do to get the solution I want with NDSolve? The reason I'm asking is because I want to solve a modified version of the equation, which DSolve can't solve but NDSolve should be able to, except I get the same problem.

Answer 1

Differentiating the ODE yields

eq = A'[t]^2 == 1/A[t] - 1;
ic = A[0] == 1;
D[eq, t] // Simplify
(*A'[t] (1/A[t]^2 + 2 A''[t]) == 0  *)

which has two solutions A'[t] == 0 and A''[t] == -1/2 1/A[t]^2, both valid at the initial condition A[0] == 1. It should be no surprise that since A[t] == 1 implies A'[t] == 0 in eq (exactly in floating point), the value of A[t] == 1 never changes in a numerical solution. This solution is numerically unstable, but since there is no round off error, the integration is never bumped off it. If there were round off error, the solution would drift until the other solution takes over.

One workaround is to use the second-order factor for the ODE, using eq to solve for A'[0]:

NDSolve[{{A''[t] == -(1/2) (1 + A'[t]^2)^2,
  A[0] == 1, A'[0] == 0},
 A, {t, 0, 2};

Another workaround is to adapt @Mariusz's idea from the comments to a value t sufficiently close to the IC but otherwise arbitrary, since this problem has an asymptotic solution at the IC. One can inspect the terms of the asymptotic expansion to see when machine precision is achieved (assuming that terms continue to decrease). In the OP's problem, order 10 and above is a high enough approximation to start integration at t == 0.1.

AsymptoticDSolveValue[{eq, ic, A''[0] == -1/2}, A, {t, 0, 12}] // 
 Apply@List
% /. t -> 0.1
(*
{1, -(t^2/4), -(t^4/48), -((11 t^6)/2880), -((73 t^8)/80640),
 -((887 t^10)/3628800), -((136883 t^12)/1916006400)}

{1, -0.0025, -2.08333*10^-6, -3.81944*10^-9,
 -9.05258*10^-12, -2.44433*10^-14, -7.14418*10^-17}
*)

Since eq is quadratic in A'[t], there are two solutions at every initial condition. We want the one with a negative square root, which is the first solution returned by Solve[]:

NDSolve[
  {First@Solve[eq, A'[t]] /. Rule -> Equal,
   A[t] == 
     AsymptoticDSolveValue[{eq, ic, A''[0] == -1/2}, 
      A, {t, 0, 12}] /. t -> 0.1`,
   (* Optional: stop integration when A'[t] becomes complex *)
   WhenEvent[! Developer`RealQ[A'[t]], "StopIntegration"]},
  A, {t, 0, 2}] /.
 (* Optional: The values have a zero imaginary part 
    but are not Real. This converts the numbers to Real.
    Requires the integration to have been stopped when
    A'[t] is complex. *)
 a_ /; ArrayQ[a, 1, Developer`MachineComplexQ] :> Re[a]

For those who inspect every detail, there is a little round off error at the end of the last step that leads to an imaginary part approximately equal to $3\times10^{-15}$, which is lopped off by Re[]. One can use the WhenEvent option "LocationMethod" -> "StepBegin", to avoid it, if it is bothersome.

Answer 2

Since Reduce[eq, A'[t], Reals] says 0 < A[t] <= 1 && ..... linearise expected solution Asol at t==0 with $MachineEpsilon. Yields a fast and quite good result with only small errors at tmax==Pi/2. (Even better at higher workingprecision)

Asol = A /. 
  First@NDSolve[{eq, A[0 + $MachineEpsilon] == 1 - $MachineEpsilon}, 
    A, {t, 0 - $MachineEpsilon, 2}, Method -> "StiffnessSwitching"]

{Asol[0], Asol[Pi/2]}

(*   {1., 0.00434685 + 0. I}   *)

Plot[Asol[t], {t, 0, 2}, PlotRange -> All]