Incorrect results when solving an ODE

Question

The solving of an ODE problem in version 13 on Windows 10

sol = DSolve[{x'[t] == Sqrt[1 - 2 x[t]^2]^2 + x[t]*Sqrt[1 - 2 x[t]^2],
x[0] == 0}, x, t]

{{x -> Function[{t}, ( 2 (Sqrt[3] Tanh[(Sqrt[3] t)/2] + Tanh[(Sqrt[3] t)/2]^2))/( 3 + 2 Sqrt[3] Tanh[(Sqrt[3] t)/2] + 3 Tanh[(Sqrt[3] t)/2]^2)]}, {x -> Function[{t}, (2 (Sqrt[3] Tanh[1/2 (Sqrt[3] t - 2 ArcTanh[Sqrt[3]])] + Tanh[1/2 (Sqrt[3] t - 2 ArcTanh[Sqrt[3]])]^2))/(3 + 2 Sqrt[3] Tanh[1/2 (Sqrt[3] t - 2 ArcTanh[Sqrt[3]])] + 3 Tanh[1/2 (Sqrt[3] t - 2 ArcTanh[Sqrt[3]])]^2)]}}

brings two solutions and this contradicts the existence and uniqueness of the solution of the Cauchy problem. The above is not a usual technical bug in Mathematica. A bunch of problems related to that ODE problem appears. First, the results of

Plot[(x'[t] - Sqrt[1 - 2 x[t]^2]^2 - x[t]*Sqrt[1 - 2 x[t]^2]) /. 
sol[[1]], {t, 0, 5}, PlotRange -> All] 

and

Plot[(x'[t] - Sqrt[1 - 2 x[t]^2]^2 - x[t]*Sqrt[1 - 2 x[t]^2]) /. sol[[2]], {t, 0, 5}]

demonstrate that sol[[1]] is not a solution of that ODE at all and sol[[2]] is its solution only on the interval $[0,t_0]$ with $t_0\approx 1.04$.

Second, the command

FullSimplify[{x'[t] - Sqrt[1 - 2 x[t]^2]^2 - x[t]*Sqrt[1 - 2 x[t]^2], 
x[0] == 0} /. sol[[2]], Assumptions -> t >= 0 && t <= 1]

is running without any response on my comp during more than two hours. Likely an infinite loop is accomplished here as well as in many, many other commands of Mathematica since the resourses of my comp are not exhausted (40 - 50 % of its processor and 200 - 350 MB of its ROM).

Third, when numerically solving the problem under consideration by

soln = NDSolve[{x'[t] == Sqrt[1 - 2 x[t]^2]^2 + x[t]*Sqrt[1 - 2 x[t]^2], x[0] == 0}, 
x, {t, 0, 5}, MaxSteps -> 10^6];
Plot[Evaluate[x[t] /. soln], {t, 0, 2}, PlotRange -> All]

, one obtains the solution only on the interval $[0,t_0]$ with $t_0\approx 1.04$.

For comparison, Maple 2022 results in x(t) = RootOf(t + Int(-1/(1 - 2*_a^2 + _a*sqrt(-2*_a^2 + 1)), _a = 0 .. _Z)) for the symbolical solution and the same result as Mathematica for the numeric solution.

I think the above problems are caused by

x[1.04] /. sol[[2]]

0.707105 - 9.74433*10^-20 I

so DSolve and NDSove are not able to go over the branch point of Sqrt[1 - 2 x[t]^2] at $t_0\approx 1.04$.

The questions arise: is there a workaround? how to continue the solution for $t \ge t_o$?

Answer 1

Edit2

Defining a piecewise function according to the result of ndsol (which is equal to analytical solution xsol = x /. DSolve[{eq, x[0] == 0}, x, t] below 1.04) proofs this to be a valid solution.

xn = Function[t, 
  Evaluate[Piecewise[{{xsol[[1]][t] // FullSimplify, 
      0 <= t <= Log[1 + Sqrt[3] + Sqrt[6 + 3 Sqrt[3]]]/Sqrt[3]}}, 
    1/Sqrt[2]]]]

Plot[Evaluate[xn[t]], {t, 0, 3}, PlotRange -> All]

Plot[Evaluate[(Subtract @@ eq) /. x -> xn], {t, 0, 3}, 
 PlotRange -> All, PlotStyle -> Thick]

I'm working with the older version "8.0 for Microsoft Windows (32-bit) (December 9, 2010)"

Another way: Splitting x[t] explicitly in real and imaginary part by x -> (a[#] + I b[#] &) Get same results, See below.

You get at least a numerical solution beyond t == Log[1 + Sqrt[3] + Sqrt[6 + 3 Sqrt[3]]]/Sqrt[3] where x[t] == 1/Sqrt[2] if you NDSolve the complexExpanded equation.

eq = x'[t] == Sqrt[1 - 2 x[t]^2]^2 + x[t]*Sqrt[1 - 2 x[t]^2];

ceeq = ComplexExpand[eq, TargetFunctions -> {Re, Im}] // Simplify

(*   2 x[t]^2 + Derivative[1][x][t] == 
 1 + (Cos[1/2 ArcTan[1 - 2 x[t]^2, 0]] + 
     I Sin[1/2 ArcTan[1 - 2 x[t]^2, 0]]) x[t] ((1 - 2 x[t]^2)^2)^(1/4)   *)

xsol = x /. 
  Flatten@NDSolve[{ceeq, x[0] == 0}, x, {t, 0, 5}, MaxSteps -> 10^6]

(*   NDSolve::mxst: Maximum number of 1000000 steps reached at the point t == 3.014453178076721`. >>   *)

Plot[{Re@xsol[t], Im@xsol[t]}, {t, 0, 3}, PlotRange -> All, 
 PlotStyle -> {Blue, {Thick, Red}}]

You see, beyond the critical point, solutions remains constant.

Test with leftsided equation confirms this solution.

Plot[Evaluate[Through[{Re, Im}[Subtract @@ eq /. x -> xsol]]], {t, 0, 
  3}, PlotRange -> All, PlotStyle -> {Blue, Red}]

Edit

eq = x'[t] == Sqrt[1 - 2 x[t]^2]^2 + x[t]*Sqrt[1 - 2 x[t]^2]

ff2 = (Subtract @@ eq) /. x -> (a[#] + I b[#] &)

(*   -1 - Sqrt[1 - 2 (a[t] + I b[t])^2] (a[t] + I b[t]) + 
 2 (a[t] + I b[t])^2 + Derivative[1][a][t] + I Derivative[1][b][t]   *)

cere = ComplexExpand[Re@ff2, TargetFunctions -> {Re, Im}] // Simplify

(*   -1 + 2 a[t]^2 - 2 b[t]^2 - 
 a[t] (16 a[t]^2 b[t]^2 + (1 - 2 a[t]^2 + 2 b[t]^2)^2)^(1/4)
   Cos[1/2 ArcTan[1 - 2 a[t]^2 + 2 b[t]^2, -4 a[t] b[t]]] + 
 b[t] (16 a[t]^2 b[t]^2 + (1 - 2 a[t]^2 + 2 b[t]^2)^2)^(1/4)
   Sin[1/2 ArcTan[1 - 2 a[t]^2 + 2 b[t]^2, -4 a[t] b[t]]] + 
 Derivative[1][a][t]   *)

ceim = ComplexExpand[Im@ff2, TargetFunctions -> {Re, Im}] // Simplify

(*   4 a[t] b[t] - 
 b[t] (16 a[t]^2 b[t]^2 + (1 - 2 a[t]^2 + 2 b[t]^2)^2)^(1/4)
   Cos[1/2 ArcTan[1 - 2 a[t]^2 + 2 b[t]^2, -4 a[t] b[t]]] - 
 a[t] (16 a[t]^2 b[t]^2 + (1 - 2 a[t]^2 + 2 b[t]^2)^2)^(1/4)
   Sin[1/2 ArcTan[1 - 2 a[t]^2 + 2 b[t]^2, -4 a[t] b[t]]] + 
 Derivative[1][b][t]   *)

ndsol = NDSolve[{cere == 0, ceim == 0, a[0] == 0, b[0] == 0}, {a, 
   b}, {t, 0, 3}, MaxSteps -> 10^6]

(*   NDSolve::mxst: Maximum number of 1000000 steps reached at the point t == 2.4880047588617566`. >>  *)

Plot[Evaluate[{a[t], b[t]} /. ndsol], {t, 0, 2.48}, PlotRange -> 1, 
 GridLines -> {{Log[1 + Sqrt[3] + Sqrt[6 + 3 Sqrt[3]]]/Sqrt[3]}, {1/
     Sqrt[2]}}, PlotStyle -> Thick]

Plot[Evaluate[Re[ff2 /. ndsol]], {t, 0, 2.48}]

Plot[Evaluate[Im[ff2 /. ndsol]], {t, 0, 2.48}]