3

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?

sara nj
  • 311
  • 1
  • 6

2 Answers2

5

Modifying my answer here, this gives a solution:

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

 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.


Update

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 // 
   First;
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[{
   ode,
   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}]

ClearAll[cp];
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]},
 Show[
  cp[v],
  Table[
   pf1[y1, p1, v] // Quiet;
   ParametricPlot3D[
     {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[
       ReplacePart[
        RotateLeft[Join @@ Table[{.0, .0, .03, .0, .0}, {5}], 
         3 + 2 Round[4 Abs[p1] + 1/4]], {1 -> 0., -1 -> 0.}]], 
      Arrow[Tube@p]},
   {p1, -1.875, 1.875, 1/4}],
  Graphics3D[{
    {cf[6],
     Scale[
      Tube[{{-Pi/2, 0, 0}, {Pi, 0, 0}}], {1, 1/2, 1/4}, {0, 0, 0}]},
    {cf[7],
     Table[
      Scale[
       Tube[{{y, 0, -4}, {y, 0, 4}}], {2/(3 Pi), 1/2, 1}, {y, 0, 0}],
      {y, 
       y /. Solve[-Pi/2 <= y <= Pi && 
          Coefficient[manifold /. Equal -> Subtract, q] == 0, y]}
      ]
     }
    }],
  Axes -> Automatic, Lighting -> "Neutral"
  ]
 ]
Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • 3
    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. – Bob Hanlon May 30 '20 at 16:10
  • 1
    @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[]. – Michael E2 Jun 03 '20 at 00:20
0

This is an extended comment rather than an answer.

Clear["Global`*"]

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}]

sol[-1][x]

enter image description here

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

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

RootApproximant[%/Pi]*Pi

(* {Pi/2, Pi/2, Pi/2, Pi/2, Pi/2} *)
Bob Hanlon
  • 157,611
  • 7
  • 77
  • 198
  • Isn't obvious that $f = \pi/2$ (constant function) is a solution? – Michael E2 May 30 '20 at 17:44
  • @MichaelE2 - I am trying to elicit a consistent set of boundary conditions for the problem. – Bob Hanlon May 30 '20 at 17:48
  • f[x] could be anything so it doesn't need to have a period of $2 \pi$ – sara nj May 30 '20 at 17:56
  • I don't know how to choose a consistent boundary condition! – sara nj May 30 '20 at 17:57
  • 1
    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. – Michael E2 May 30 '20 at 18:02
  • But a constant solution is very obvious! – sara nj May 30 '20 at 18:06