8
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ NDSolve give you only singular solution. Try with A[Pi/2]=0.0001. Dosen't work with: A[Pi/2]=0 very strange !. $\endgroup$ Aug 27, 2022 at 15:32
  • $\begingroup$ Your solution works, but only because thanks to DSolve we can know some values at other points. For the modified equation, I only know that the initial value A[0] is 1 and can't plug in other values to try and "trick" NDSolve into giving me the correct solution. $\endgroup$ Aug 27, 2022 at 15:59
  • $\begingroup$ You can differentiate the ODE a number of times and reduce it to eliminate divide-by-zero problems (connected with the singular solution) at the IC; and use the differentiated equations to extend the IC to the needed extra values $A^{(k)}(0) = a^{(k)}_0$. For instance, applying NDSolve to {A''[t] == -(1/2) (1 + A'[t]^2)^2, A[0] == 1, A'[0] == 0} yields the desired solution. $\endgroup$
    – Michael E2
    Aug 27, 2022 at 18:09

2 Answers 2

4
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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]
$\endgroup$
2
  • $\begingroup$ Note that A'[0] == 0 if A[0] == 1 so instead of linearizing, since the linear term is zero, you need to go out to the quadratic term or higher. I went higher to move the IC further away from t == 0, so as to avoid loss of precision due to round-off. And AsymptoticDSolveValue[] makes this easy. $\endgroup$
    – Michael E2
    Aug 28, 2022 at 13:31
  • $\begingroup$ @MichaelE2 , you are right. Shouldn't call it linearization. It is just going smallest step away from zero and supposing, solution A to be the same step smaller like y==1-x. $\endgroup$
    – Akku14
    Aug 28, 2022 at 14:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.