2
$\begingroup$

Introduction and problem

I tried to solve, symbolically and numerically, the following ODE with Mathematica: $y’(t) = \sin(t) y’’(t)$, and I got strange results at different levels (either inconsistent, or errors).

I point out first that I approached this problem having already found a particular solution by hand: $y(t) = \log(\cos(\frac{x}{2}))$.

Inconsistent results between versions with DSolve

I first tried to solve this equation with Mathematica 12:

$Version
(* 12.3.1 for Microsoft Windows (64-bit) (June 24, 2021) *)

DSolve[y'[t] == Sin[t] y''[t], y[t], t]
(* {{y[t] -> C[2] - 2 C[1] Log[Cos[t/2]]}} *)

I thus obtain a closed-form solution, more general but compatible with my particular solution (with $c_1 = - \frac{1}{2}$ and $c_2 = 0$). Everything looks good so far. But now let us try to solve the very same equation, on the latest version Mathematica 14:

$Version
(* 14.0.0 for Linux x86 (64-bit) (December 13, 2023) *)

DSolve[y'[t] == Sin[t] y''[t], y[t], t]
(* {{y[t] -> C[2] + C[1] Csc[t] Log[Sec[t/2]^2] Sqrt[Sin[t]^2]}} *)

Well, we get a slightly different solution… One might think it’s another form of the same expression, or an even more general one, but it’s not.

Firstly, it is defined in more places than the previous form if we assume that the function is real (i.e. in places where $\cos(\frac{t}{2})$ was negative, $\log(\cos(\frac{t}{2}))$ was not defined).

And even if we can see that for certain values the two forms simplify and are equal, a simple study of the signs shows that they are sometimes of opposite signs for certain values, regardless of the constant values $c_1$, $c_2$, i.e. when $4 \pi n - \pi < t < 4 \pi n$ for all $n \in \mathbb{Z}$. The two forms can also be plotted to see that: Plot of the two forms between −5 and 15

I’ve checked, the two forms returned seem to be valid solutions; but for me there’s a problem, because I expect DSolve to return all possible solutions, and not to have different and incompatible solutions depending on the version of Mathematica used.

Impossible to solve the equation numerically with NDSolve

Finally, I also tried to solve the same equation numerically between $-1$ and $1$ (an interval where it is numerically stable), but this time I get an error:

NDSolve[{y'[t] == Sin[t] y''[t], y[0] == 0, y'[0] == 0}, y, {t, -1, 1}]

(* Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. *)
(* NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0. *)
(* NDSolve[{Derivative[1][y][t] == Sin[t] (y^\[Prime]\[Prime])[t], y[0] == 0, Derivative[1][y][0] == 0}, y, {t, -1, 1}] *)

And this time, all versions of Mathematica give the same result 😅.

Conclusion

I think I’ve found two cases where Mathematica’s outputs are incomplete or missing, and would therefore be bugs. But before reporting bug reports to Wolfram, I wanted to discuss these results to better understand what I’m getting, and above all to check that it’s not me who’s done something wrong.

$\endgroup$
5
  • 1
    $\begingroup$ Re "if we assume that the function is real": Mathematica is assuming the functions and constants C[1], C[2] are complex. Further, equivalence need only exists over intervals that do not contain singular points ($t=\pi n$, $n \in \Bbb Z$). Over different intervals, the two forms might be related by different constants. Your NDSolve code has the initial condition at a singular point, which is why it fails. $\endgroup$
    – Goofy
    Commented Mar 17 at 14:01
  • $\begingroup$ In your solution with NDSolve[] the integration limits are the same, thus the code fails. $\endgroup$
    – Hans Olo
    Commented Mar 17 at 15:10
  • $\begingroup$ But I'm not sure if this is the only issue at the moment (writing on mobile)... $\endgroup$
    – Hans Olo
    Commented Mar 17 at 15:10
  • 1
    $\begingroup$ @HansOlo While the OP clearly meant {t, -1, 1}, it's not an error for the end points to be the same. Try out this: sol = NDSolve[{y'[x] == y[x], y[0] == 1}, y, {x, 1, 1}]; y[1] /. sol. $\endgroup$
    – Goofy
    Commented Mar 17 at 15:30
  • $\begingroup$ @HansOlo it was a typo, Goofy is right, fixed! $\endgroup$
    – mlpo
    Commented Mar 17 at 15:34

1 Answer 1

2
$\begingroup$

To illustrate the first part of my comment:

Re "if we assume that the function is real": Mathematica is assuming the functions and constants C[1], C[2] are complex. Further, equivalence need only exists over intervals that do not contain singular points ($t=\pi n$, $n\in\Bbb Z$). Over different intervals, the two forms might be related by different constants. Your NDSolve code has the initial condition at a singular point, which is why it fails.

It is important to note that the constants of integration of the two solutions are independent. Thus I will change the second solution's constants to C[3], C[4].

ysol12 = {{y[t] -> C[2] - 2  C[1]  Log[Cos[t/2]]}} // First // 
  DSolve`DSolveToPureFunction
ysol14 = {{y[t] -> 
       C[2] + C[1]  Csc[t]  Log[Sec[t/2]^2]  Sqrt[Sin[t]^2]}} /. 
    C -> (C[# + 2] &) // First // DSolve`DSolveToPureFunction
(*
{y -> Function[{t}, C[2] - 2  C[1]  Log[Cos[t/2]]]}

{y -> Function[{t}, C[4] + C[3] Csc[t] Log[Sec[t/2]^2] Sqrt[Sin[t]^2]]}
*)

Now we can solve for the relationship between the pairs of constants that make the solutions equivalent over various intervals:

(* 0 < t < Pi *)
{y[1], y[2]} /. {ysol12, ysol14} // Apply@Equal // 
  Solve // FullSimplify
(* {{C[3] -> C[1], C[4] -> C[2]}} *)

(* Pi < t < 2Pi *)
{y[4], y[5]} /. {ysol12, ysol14} // Apply@Equal // 
  Solve // FullSimplify
(* {{C[3] -> -C[1], C[4] -> -2  I  \[Pi]  C[1] + C[2]}} *)

(* 2Pi < t < 3Pi *)
{y[7], y[9]} /. {ysol12, ysol14} // Apply@Equal // 
  Solve // FullSimplify
(* {{C[3] -> C[1], C[4] -> -2  I  \[Pi]  C[1] + C[2]}} *)

The problem with NDSolve is simple mistake. It's not clear whether further comment is needed. One might point out that because there is some error at each step (usually small), the solver solves a new IVP at each step. When the solution approaches a singular point, then the solver might stop or it might step over the singularity. When it passes a singular point, then the error might be quite large (and undetected). Keep in mind that there isn't necessarily a unique solution at a singular point, and the solver might step onto any of one of the solutions. Further, which solution is followed may depend on the integration method, the PrecisionGoal, the WorkingPrecision, etc. For instance, in the OP's ODE, y'[t] is forced to approach $0$ as t approaches 0. Any error in y'[t] can represent an arbitrary loss of precision, since y'[t] is approximating zero at the singular point. And that's not the fault of NDSolve.

$\endgroup$
3
  • $\begingroup$ Thanks for the comprehensive answer! For symbolic forms, I didn’t have in mind that “equivalence need only exists over intervals that do not contain singular points” (and although this may seem logical, it doesn’t seem to me to be very well documented in the DSolve page 😅). For NDSolve, I indeed hadn’t taken the singular point into account. Here, by giving initial conditions at $t = \frac{1}{2}$, instead of $t = 0$, to avoid the singular point, I get a valid numeric solution. Are there any other simple ways of preventing singular points from getting in the way of the solver? $\endgroup$
    – mlpo
    Commented Mar 18 at 15:44
  • 1
    $\begingroup$ @mlpo It's a little glib, but one might say that singular points are why there are graduate schools. Intro undergrad DE courses, if they go beyond saying that at singular points you lose the guarantee of existence and uniqueness of solutions, discuss the Frobenius series method. And they might not point out that these solutions might not be defined on both sides of the singular point. For instance, $y =\cases{ a x^2 & x<0\cr 0 & x=0\cr b x^2 & x>0\cr}$ is the most general solution to $x y' = 2y$ for any $a$ and $b$. There is no theory that says we must have $a=b$... $\endgroup$
    – Goofy
    Commented Mar 18 at 16:16
  • $\begingroup$ ...When the solver reaches a singular point, there may be infinitely many solutions or none. If there are many, then any one it computes is correct. I called it an "error," but there is no theoretical exact continuation to compare it to. -- I suppose you could require an analytic solution (that is, a function that can be represented by a power series). Then you would have uniqueness, if a solution existed (it might not). $\endgroup$
    – Goofy
    Commented Mar 18 at 16:16

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.