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.
NDSolve
give you only singular solution. Try withA[Pi/2]=0.0001
. Dosen't work with:A[Pi/2]=0
very strange !. $\endgroup$NDSolve
to{A''[t] == -(1/2) (1 + A'[t]^2)^2, A[0] == 1, A'[0] == 0}
yields the desired solution. $\endgroup$