I have found the solution of an ODE in Mathematica by the below script:

DSolve[2 D[f[x], x]^2 (-v Cos[2 f[x]] + Sin[2 f[x]]) - 
   D[D[f[x], x], x] (Cos[2 f[x]] + v Sin[2 f[x]]) == 0, f[x], x]

I'm wondering if it is possible to constrain Mathematica to only print solutions that are real for real values of $x$.Could someone please help?

  • $\begingroup$ What are the boundary conditions and what if any constraints exist on v? $\endgroup$
    – Bob Hanlon
    Commented May 30, 2020 at 15:27
  • $\begingroup$ Without another boundary condition there is still a generated parameter. And the additional boundary condition must be consistent. $\endgroup$
    – Bob Hanlon
    Commented May 30, 2020 at 16:20

2 Answers 2


Modifying my answer here, this gives a solution:

Internal`InheritedBlock[{Solve},   (* hijack Solve to give only Real sols *)
 Solve[eq_, v_, opts___] /; ! TrueQ[$in] := 
  Block[{$in = True, $res1, $res2}, 
   Solve[eq, v, Reals, Method -> Reduce, opts]];

 DSolve[2 D[f[x], x]^2 (-v Cos[2 f[x]] + Sin[2 f[x]]) - 
    D[D[f[x], x], x] (Cos[2 f[x]] + v Sin[2 f[x]]) == 0, f[x], x]


The result is sufficiently complicated that I will let the OP verify it is what is desired.


This is a more elaborate look at the ODE than is necessary, mainly because ODEs are fun and this turns out to be a nice example. I've written before on the difficulties posed by a distinct component of the solution space that is not the limit/boundary of the general solution in the first link below; the second link is only somewhat related, primarily by its containing a similar form of analysis as I will present below. This look at the ODE will also give some insight into feasible BVPs, which alluded to in a comment.

The "contact manifold" for a second-order autonomous ODE can be thought of as a surface in space with the coordinates $y=f(x)$, $p=f'(x)$, $q=f''(x)$:

ode = 2 D[f[x], x]^2 (-v Cos[2 f[x]] + Sin[2 f[x]]) - 
    D[D[f[x], x], x] (Cos[2 f[x]] + v Sin[2 f[x]]) == 0;
manifold = ode /. {f''[x] -> q, f'[x] -> p, f[x] -> y}

(*  2 p^2 (-v Cos[2 y] + Sin[2 y]) - q (Cos[2 y] + v Sin[2 y]) == 0  *)

The ODE/manifold has a cute representation, of which I cannot figure out how to take advantage: $$\pmatrix{2p^2 & q} \pmatrix{ \cos 2 y & -\sin 2 y \\ \sin 2 y & \phantom{-}\cos 2 y \\} \pmatrix{v \\ 1\\} = 0 \,. $$

There is a natural direction field on this surface arising from the relation $q = d^2y/dx^2 = p \, dp/dy$. We get a system of "contact planes" $q\, dy = p\, dp$; more explicitly, at each point $X_0=(y_0,p_0,q_0)$ with $(p_0,q_0)\ne(0,0)$, we have the plane $q_0 (y-y_0) = p_0 (p-p_0)$. Where the planes intersect the manifold transversely at $X_0$ (i.e., the contact plane and tangent plane intersect in a line), they define a tangent direction. Thus the system defines a direction field on the contact manifold. The trajectories of the solutions to the ODE will be tangent to this "contact field." This makes a very helpful way to visualize the solution system of the ODE. For the OP's example, it will be worth noting that the contact planes along the $y$ axis, where $p=q=0$, are undefined. It turns out that the contact manifold in this case contains the $y$ axis, so it forms a singular locus where the contact field is undefined. One will note further that there are singular points of the ODE where the coefficient of $q$ is zero. They are at $$y_s = \tan ^{-1}\left(v\pm\sqrt{v^2+1}\right) + \pi n = \tan ^{-1}\left(v+\sqrt{v^2+1}\right) + {\pi\over2}\, n, \quad n \in {\Bbb Z}\,. \tag{1}$$ For these values of $y = y_s$, the only IVPs with solutions to the given form have $p=0$.

Three views of the contact manifold

Fig. 1. Three views of the contact manifold (v = 1). The flow of the contact field is illustrated on one of the sheets. The red line is the singular locus (the $y$-axis) and corresponds to the constant solutions of the ODE. Note that the flow as it approaches the red line becomes parallel to it.

One can see that the non-constant solutions are monotonic, because the only solutions that can have a zero derivative $p$ are the constant solutions:

Reduce[{manifold, p == 0, y ∈ Reals}]
  ((y ∈ Reals && q == 0) ||                       <-- y-axis
   (v ∈ Reals && C[1] ∈ Integers &&
     (y == ArcTan[v - Sqrt[1 + v^2]] + π C[1] ||  <-- vert lines
      y == ArcTan[v + Sqrt[1 + v^2]] + π C[1]))       at sing pts
   ) && p == 0

It is clear that a necessary condition for a BVP $f(x_1)=y_1,\ f(x_2)=y_2$ to have a solution is that $y1, y2$ must lie between two consecutive singular points $y_s$ in (1). It is also clear that a solution, if it exists, is unique and that if $y_1 = y_2$, then the constant solution is the only solution. Since for any $y \ne y_s$, the rate $p$ can be as close to zero or as large in magnitude as we like, it seems likely that

Every BVP of the OP's ODE has a (unique) solution.

Solving any IVP

The OP's ODE can be solved by hand (I'll use Integrate to help because this is mma.SE, not math.SE) for any IVP $f(x_0)=y_0$, $f(x_0)=p_0$ as long as $y_0 \ne y_s$ or $p_0 \ne 0$. It will be seen below that these two conditions come from the solution comprising two components.

As we did for the contact field, I will substitute $q = p \, dp/dy$, which yields

Factor /@ (manifold /. q -> p p')
(* -p (2 (v Cos[2y] - Sin[2y]) p + (Cos[2y] + v Sin[2y]) p') == 0  *)

This splits into two components, $$ p=0, \quad (v \sin 2 y + \cos 2 y)\, {dp \over dy} + 2 (v \cos 2 y - \sin 2y)\,p = 0 \,.$$ The first $p = dy/dx = 0$ yields the constant functions and the second yields the general solution returned by DSolve. The second one is separable and so can be integrated twice. We'll use an initial condition $(x_0,y_0,p_0)$ (i.e., $f(x_0)=y_0$, $f'(x_0)=p_0$). Recall that if $g(p) \,dp = h(y)\,dy$, then the solution to the IVP satisfies $$\int_{p_0}^p g(p) \; dp = \int_{y_0}^y h(y) \; dy \,.$$ The first integration yields (for $p_0\ne0$) $$\log p - \log p_0 = \log(v \sin 2y_0 + \cos 2y_0) - \log(v \sin 2y + \cos 2y) \,,$$ which is equivalent to $$ {1 \over p_0}\,{dy \over dx} = { v \sin 2y_0 + \cos 2y_0 \over v \sin 2y + \cos 2y }$$ and easily integrated again.

The code to get this is the following:

rhs1 = Integrate[-((2 (v Cos[2 y] - Sin[2 y]))/(
         Cos[2 y] + v Sin[2 y])), y] /. {{y -> y0}, {y -> y}} // 
      Differences // First // Exp // Simplify;
lhs1 = Integrate[1/p, p] /. {{p -> p0}, {p -> p}} // Differences // 
     First // Exp // Simplify;
p /. First@Solve[lhs1 == rhs1, p]
(*  (p0 (Cos[2 y0] + v Sin[2 y0]))/(Cos[2 y] + v Sin[2 y])  *)

The second integration yields an implicit solution:

rhs2 = Integrate[(p0 (Cos[2 y0] + v Sin[2 y0])/(
          Cos[2 y] + v Sin[2 y]))^-1, y] /. {{y -> y0}, {y -> y}} // 
     Differences // First // Simplify;
lhs2 = Integrate[1, x] /. {{x -> x0}, {x -> x}} // Differences // 
lhs2 == rhs2
  x - x0 ==
    (Sin[y - y0] (Cos[y + y0] + v Sin[y + y0]))/
     (p0 (Cos[2 y0] + v Sin[2 y0]))

The solution of this (not particularly easy to obtain with Mathematica) gives the DSolve solution. The implicit solution can be verified with the following:

fpsol = First@Solve[D[lhs2 == rhs2 /. y -> f[x], x], f'[x]];
fppsol = First@Solve[D[lhs2 == rhs2 /. y -> f[x], x, x], f''[x]];
ode /. fppsol /. fpsol // Simplify

(*  True  *)

Code dump for graphics

 cf = ColorDataFunction["OkabeIto", "Indexed", {1, 8, 1},
     {RGBColor[0.902, 0.624, 0], RGBColor[0.337, 0.706, 0.914], 
      RGBColor[0, 0.62, 0.451], RGBColor[0.941, 0.894, 0.259], 
      RGBColor[0, 0.447, 0.698], RGBColor[0.835, 0.369, 0], 
      RGBColor[0.8, 0.475, 0.655]}[[Mod[Floor[#1], 8, 1]]] &];

pf1 = ParametricNDSolveValue[{
   f[0] == y1(*Pi+ArcTan[(-1-(v Sqrt[1+v^2])/Abs[v])/v]*),
   f'[0] == p1,
   WhenEvent[Abs[f''[x]] > 4, "StopIntegration"]
   }, f, {x, -50, 50}, {y1, p1, v},
  "ExtrapolationHandler" -> {Indeterminate &, 
    "WarningMessage" -> False}]

mem : cp[v0_] := mem = ContourPlot3D[
    Evaluate[manifold /. v -> v0],
    {y, -Pi/2, Pi}, {p, -2, 2}, {q, -4, 4},
    Axes -> Automatic, AxesLabel -> Automatic, Mesh -> None, 
    ContourStyle -> {Opacity[0.8], FaceForm[cf[4]]}

Block[{v = 1/2, y1 = Pi + ArcTan[(-1 - (v Sqrt[1 + v^2])/Abs[v])/v]},
   pf1[y1, p1, v] // Quiet;
     {pf1[y1, p1, v][x], pf1[y1, p1, v]'[x], pf1[y1, p1, v]''[x]},
     Evaluate@Flatten@{x, pf1[y1, p1, v]["Domain"]},
     PlotStyle -> cf[5]] /. 
    Line[p_] :> {Arrowheads[
        RotateLeft[Join @@ Table[{.0, .0, .03, .0, .0}, {5}], 
         3 + 2 Round[4 Abs[p1] + 1/4]], {1 -> 0., -1 -> 0.}]], 
   {p1, -1.875, 1.875, 1/4}],
      Tube[{{-Pi/2, 0, 0}, {Pi, 0, 0}}], {1, 1/2, 1/4}, {0, 0, 0}]},
       Tube[{{y, 0, -4}, {y, 0, 4}}], {2/(3 Pi), 1/2, 1}, {y, 0, 0}],
       y /. Solve[-Pi/2 <= y <= Pi && 
          Coefficient[manifold /. Equal -> Subtract, q] == 0, y]}
  Axes -> Automatic, Lighting -> "Neutral"
  • 3
    $\begingroup$ This does not work with her boundary condition (given in a comment) of f[0] == Pi/2. I suspect that restricting the domain of Solve to Reals precludes necessary intermediate results containing imaginary components. $\endgroup$
    – Bob Hanlon
    Commented May 30, 2020 at 16:10
  • 1
    $\begingroup$ @BobHanlon This happens a lot, esp. with solutions expressed with inv. trig. For example, you cannot get Pi/2 from ArcTan[m/n]; however, you can if you replace it by ArcCot[n/m]. That approach will work in this case. Sometimes taking the Tan[] can solve it, if you can isolate the ArcTan[] term. A similar problem sometimes arises with Log[m/n] when n -> 0; in that case, one might apply Exp[]. $\endgroup$
    – Michael E2
    Commented Jun 3, 2020 at 0:20

This is an extended comment rather than an answer.


f[0] cannot be multivalued so Mod[f[0], 2Pi]doesn't add much. If instead you believe f[x] has a period of 2 Pi perhaps you mean f[2 n Pi] == Pi/2 where Element[n, Integers]. However, then f[x] is the constant Pi/2

sol[v_?NumericQ] := 
 NDSolveValue[{2 D[f[x], x]^2 (-v Cos[2 f[x]] + Sin[2 f[x]]) - 
     D[D[f[x], x], x] (Cos[2 f[x]] + v Sin[2 f[x]]) == 0, f[0] == Pi/2, 
   f[2 Pi] == Pi/2}, f, {x, 0, 2 Pi}]


enter image description here

Table[sol[v][RandomReal[2 Pi]], {v, -5, -1}]

(* {1.5708, 1.5708, 1.5708, 1.5708, 1.5708} *)


(* {Pi/2, Pi/2, Pi/2, Pi/2, Pi/2} *)
  • $\begingroup$ Isn't obvious that $f = \pi/2$ (constant function) is a solution? $\endgroup$
    – Michael E2
    Commented May 30, 2020 at 17:44
  • $\begingroup$ @MichaelE2 - I am trying to elicit a consistent set of boundary conditions for the problem. $\endgroup$
    – Bob Hanlon
    Commented May 30, 2020 at 17:48
  • $\begingroup$ f[x] could be anything so it doesn't need to have a period of $2 \pi$ $\endgroup$
    – sara nj
    Commented May 30, 2020 at 17:56
  • $\begingroup$ I don't know how to choose a consistent boundary condition! $\endgroup$
    – sara nj
    Commented May 30, 2020 at 17:57
  • 1
    $\begingroup$ Well, f = C (any constant) is a solution. My point is about the formulation of the problem. In the OP, we're looking for a generic solution and all is ok; my answer probably works. The condition added in the comments, which is an incomplete initial condition, is problematic. The solutions $f=C$ lie on a singular locus. DSolve does not usually find those. $\endgroup$
    – Michael E2
    Commented May 30, 2020 at 18:02

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.