0
$\begingroup$

I've been trying to solve a set of differential equations that involve the square of a derivative, similar to another post I made here, and trigonometric functions, as seen here \begin{equation} \frac{d\tau}{ds} = \sin\sigma \left( E \cos\sigma - p_\tau \sin \sigma \right) \end{equation} \begin{equation} \left(\frac{d\sigma}{ds}\right)^2 = \sin^2\sigma \left[ \left( E \cos\sigma - p_\tau \sin \sigma \right)^2 - M^2 \right] \end{equation}

I followed an entirely analogous process to a suggestion that was given in that post, but I was faced with two problems:

  1. Using the boundary conditions I chose, I only get a plot that ranges in $\tau$ from (approximately) $[-1,5]$, when the expected result covers all $\tau$.

  2. I have only been able to solve the equations and get a plot for $E=0$. If $E \neq 0$, then DSolve runs for several hours on end (and is currently still doing so).

The code I implemented was the following(in the code $y$ corresponds to $\sigma$ and $t$ to $\tau$):

(*Boundary Conditions*)
E0 = 0;
M = 5;

sol = DSolve[{t'[y]^2 == (E0*Cos[y] - p*Sin[y])^2/((E0*Cos[y] - p*Sin[y])^2 - M^2), t[1] == 2}, t, y]

Table[Plot[Evaluate[t[y] /. sol], {y, -1, Pi}, PlotLabel -> p], {p, 
  5.9, 6, 0.001}]

Plot[Evaluate[t[y] /. sol /. p -> 5.942], {y, -1, Pi}, 
 FrameLabel -> {"y", "t"}, PlotStyle -> Red, Frame -> True]

The expected results are:

enter image description here

Where left corresponds to $E=0$ and right to $E^2 > M^2$.

Thank you very much for the help.

$\endgroup$
2
  • 1
    $\begingroup$ What green and blue lines mean on these pictures? $\endgroup$ Commented Sep 15, 2023 at 12:02
  • $\begingroup$ The blue lines can be ignored, they are there only as visual aid, but aren't part of the solution. The green lines correspond to the boundaries $\sigma = 0, \pi$. $\endgroup$
    – MarcelRomp
    Commented Sep 15, 2023 at 12:38

2 Answers 2

2
$\begingroup$

In this post we answer the question how to extend periodic solution t[y] up to $4 \pi$ and how to compute solution in a case of E^2>M^2. Solution in a case of E=0 is given by

E0 = 0;
M = 5;

sol = DSolve[{t'[
      y]^2 == (E0*Cos[y] - p*Sin[y])^2/((E0*Cos[y] - p*Sin[y])^2 - 
       M^2), t[1] == 2}, t, y]

Visualization

cp = Plot[Evaluate[t[y] /. sol /. p -> 5.942], {y, -Pi, Pi}, 
  FrameLabel -> {"y", "t"}, PlotStyle -> Red, Frame -> True, 
  PlotLabel -> Row[{"E = ", E0}]];

fig1 = Show[
  Show[MapAt[Translate[#, {0, 2 Pi}] &, cp, 1], cp, 
   MapAt[Translate[#, {0, -2 Pi}] &, cp, 1], PlotRange -> All], 
  Graphics[{Green, Thick, Line[{{0, -8}, {0, 13}}], 
    Line[{{Pi, -8}, {Pi, 13}}]}]]

Figure 1

Solution in a case of E0=6 is given by

E0 = 6;
M = 5;

sol1 = DSolve[{t'[
      y]^2 == (E0*Cos[y] - p*Sin[y])^2/((E0*Cos[y] - p*Sin[y])^2 - 
       M^2), t[1] == 2}, t, y]
(*{{t -> Function[{y}, (144 Cos[1] Cos[y] Sqrt[
        14 - p^2 - 36 Cos[2] + p^2 Cos[2] + 12 p Sin[2]] Sqrt[
        14 - p^2 - 36 Cos[2 y] + p^2 Cos[2 y] + 12 p Sin[2 y]] - 
       24 p Cos[y] Sin[1] Sqrt[
        14 - p^2 - 36 Cos[2] + p^2 Cos[2] + 12 p Sin[2]] Sqrt[
        14 - p^2 - 36 Cos[2 y] + p^2 Cos[2 y] + 12 p Sin[2 y]] - 
       24 p Cos[1] Sqrt[
        14 - p^2 - 36 Cos[2] + p^2 Cos[2] + 12 p Sin[2]] Sin[y] Sqrt[
        14 - p^2 - 36 Cos[2 y] + p^2 Cos[2 y] + 12 p Sin[2 y]] + 
       4 p^2 Sin[1] Sqrt[
        14 - p^2 - 36 Cos[2] + p^2 Cos[2] + 12 p Sin[2]] Sin[y] Sqrt[
        14 - p^2 - 36 Cos[2 y] + p^2 Cos[2 y] + 12 p Sin[2 y]] + 
       6 Sqrt[2]
         ArcTanh[(-11 - 
          6 I p + (-25 I + 6 p + I p^2) Tan[1])/((-6 I + 
            p) Sqrt[-11 + 12 p Tan[1] - (-25 + p^2) Tan[1]^2])] Cos[
         1] Cos[y] Sqrt[-(6 Cos[1] - p Sin[1])^2] Sqrt[
        14 - p^2 - 36 Cos[2 y] + p^2 Cos[2 y] + 12 p Sin[2 y]]
         Sqrt[-11 + 12 p Tan[1] + 25 Tan[1]^2 - p^2 Tan[1]^2] - 
       6 Sqrt[2]
         ArcTanh[(
         11 - 6 I p + 
          I (-25 + 6 I p + p^2) Tan[1])/((6 I + p) Sqrt[-11 + 
           12 p Tan[1] - (-25 + p^2) Tan[1]^2])] Cos[1] Cos[
         y] Sqrt[-(6 Cos[1] - p Sin[1])^2] Sqrt[
        14 - p^2 - 36 Cos[2 y] + p^2 Cos[2 y] + 12 p Sin[2 y]]
         Sqrt[-11 + 12 p Tan[1] + 25 Tan[1]^2 - p^2 Tan[1]^2] - 
       Sqrt[2] p ArcTanh[(-11 - 
          6 I p + (-25 I + 6 p + I p^2) Tan[1])/((-6 I + 
            p) Sqrt[-11 + 12 p Tan[1] - (-25 + p^2) Tan[1]^2])] Cos[
         1] Sqrt[-(6 Cos[1] - p Sin[1])^2] Sin[y] Sqrt[
        14 - p^2 - 36 Cos[2 y] + p^2 Cos[2 y] + 12 p Sin[2 y]]
         Sqrt[-11 + 12 p Tan[1] + 25 Tan[1]^2 - p^2 Tan[1]^2] + 
       Sqrt[2] p ArcTanh[(
         11 - 6 I p + 
          I (-25 + 6 I p + p^2) Tan[1])/((6 I + p) Sqrt[-11 + 
           12 p Tan[1] - (-25 + p^2) Tan[1]^2])] Cos[
         1] Sqrt[-(6 Cos[1] - p Sin[1])^2] Sin[y] Sqrt[
        14 - p^2 - 36 Cos[2 y] + p^2 Cos[2 y] + 12 p Sin[2 y]]
         Sqrt[-11 + 12 p Tan[1] + 25 Tan[1]^2 - p^2 Tan[1]^2] - 
       6 Sqrt[2]
         ArcTanh[(-11 - 
          6 I p + (-25 I + 6 p + I p^2) Tan[y])/((-6 I + 
            p) Sqrt[-11 + 12 p Tan[y] - (-25 + p^2) Tan[y]^2])] Cos[
         1] Cos[y] Sqrt[
        14 - p^2 - 36 Cos[2] + p^2 Cos[2] + 12 p Sin[2]]
         Sqrt[-(6 Cos[y] - p Sin[y])^2]
         Sqrt[-11 + 12 p Tan[y] + 25 Tan[y]^2 - p^2 Tan[y]^2] + 
       6 Sqrt[2]
         ArcTanh[(
         11 - 6 I p + 
          I (-25 + 6 I p + p^2) Tan[y])/((6 I + p) Sqrt[-11 + 
           12 p Tan[y] - (-25 + p^2) Tan[y]^2])] Cos[1] Cos[y] Sqrt[
        14 - p^2 - 36 Cos[2] + p^2 Cos[2] + 12 p Sin[2]]
         Sqrt[-(6 Cos[y] - p Sin[y])^2]
         Sqrt[-11 + 12 p Tan[y] + 25 Tan[y]^2 - p^2 Tan[y]^2] + 
       Sqrt[2] p ArcTanh[(-11 - 
          6 I p + (-25 I + 6 p + I p^2) Tan[y])/((-6 I + 
            p) Sqrt[-11 + 12 p Tan[y] - (-25 + p^2) Tan[y]^2])] Cos[
         y] Sin[1] Sqrt[
        14 - p^2 - 36 Cos[2] + p^2 Cos[2] + 12 p Sin[2]]
         Sqrt[-(6 Cos[y] - p Sin[y])^2]
         Sqrt[-11 + 12 p Tan[y] + 25 Tan[y]^2 - p^2 Tan[y]^2] - 
       Sqrt[2] p ArcTanh[(
         11 - 6 I p + 
          I (-25 + 6 I p + p^2) Tan[y])/((6 I + p) Sqrt[-11 + 
           12 p Tan[y] - (-25 + p^2) Tan[y]^2])] Cos[y] Sin[1] Sqrt[
        14 - p^2 - 36 Cos[2] + p^2 Cos[2] + 12 p Sin[2]]
         Sqrt[-(6 Cos[y] - p Sin[y])^2]
         Sqrt[-11 + 12 p Tan[y] + 25 Tan[y]^2 - 
         p^2 Tan[y]^2])/(2 (-6 Cos[1] + p Sin[1]) Sqrt[
       14 - p^2 - 36 Cos[2] + p^2 Cos[2] + 
        12 p Sin[2]] (-6 Cos[y] + p Sin[y]) Sqrt[
       14 - p^2 - 36 Cos[2 y] + p^2 Cos[2 y] + 
        12 p Sin[2 y]])]}, {t -> 
   Function[{y}, (-144 Cos[1] Cos[y] + 24 p Cos[y] Sin[1] + 
       24 p Cos[1] Sin[y] - 4 p^2 Sin[1] Sin[y] + 
       6 Sqrt[2]
         ArcTanh[(-11 - 
          6 I p + (-25 I + 6 p + I p^2) Tan[1])/((-6 I + 
            p) Sqrt[-11 + 12 p Tan[1] - (-25 + p^2) Tan[1]^2])] Cos[
         1] Cos[y] Sqrt[-((6 Cos[1] - p Sin[1])^2/(
         14 - p^2 - 36 Cos[2] + p^2 Cos[2] + 12 p Sin[2]))]
         Sqrt[-11 + 12 p Tan[1] + 25 Tan[1]^2 - p^2 Tan[1]^2] - 
       6 Sqrt[2]
         ArcTanh[(
         11 - 6 I p + 
          I (-25 + 6 I p + p^2) Tan[1])/((6 I + p) Sqrt[-11 + 
           12 p Tan[1] - (-25 + p^2) Tan[1]^2])] Cos[1] Cos[
         y] Sqrt[-((6 Cos[1] - p Sin[1])^2/(
         14 - p^2 - 36 Cos[2] + p^2 Cos[2] + 12 p Sin[2]))]
         Sqrt[-11 + 12 p Tan[1] + 25 Tan[1]^2 - p^2 Tan[1]^2] - 
       Sqrt[2] p ArcTanh[(-11 - 
          6 I p + (-25 I + 6 p + I p^2) Tan[1])/((-6 I + 
            p) Sqrt[-11 + 12 p Tan[1] - (-25 + p^2) Tan[1]^2])] Cos[
         1] Sqrt[-((6 Cos[1] - p Sin[1])^2/(
         14 - p^2 - 36 Cos[2] + p^2 Cos[2] + 12 p Sin[2]))]
         Sin[y] Sqrt[-11 + 12 p Tan[1] + 25 Tan[1]^2 - p^2 Tan[1]^2] +
        Sqrt[2] p ArcTanh[(
         11 - 6 I p + 
          I (-25 + 6 I p + p^2) Tan[1])/((6 I + p) Sqrt[-11 + 
           12 p Tan[1] - (-25 + p^2) Tan[1]^2])] Cos[
         1] Sqrt[-((6 Cos[1] - p Sin[1])^2/(
         14 - p^2 - 36 Cos[2] + p^2 Cos[2] + 12 p Sin[2]))]
         Sin[y] Sqrt[-11 + 12 p Tan[1] + 25 Tan[1]^2 - 
         p^2 Tan[1]^2] - 
       6 Sqrt[2]
         ArcTanh[(-11 - 
          6 I p + (-25 I + 6 p + I p^2) Tan[y])/((-6 I + 
            p) Sqrt[-11 + 12 p Tan[y] - (-25 + p^2) Tan[y]^2])] Cos[
         1] Cos[y] Sqrt[-((6 Cos[y] - p Sin[y])^2/(
         14 - p^2 - 36 Cos[2 y] + p^2 Cos[2 y] + 12 p Sin[2 y]))]
         Sqrt[-11 + 12 p Tan[y] + 25 Tan[y]^2 - p^2 Tan[y]^2] + 
       6 Sqrt[2]
         ArcTanh[(
         11 - 6 I p + 
          I (-25 + 6 I p + p^2) Tan[y])/((6 I + p) Sqrt[-11 + 
           12 p Tan[y] - (-25 + p^2) Tan[y]^2])] Cos[1] Cos[
         y] Sqrt[-((6 Cos[y] - p Sin[y])^2/(
         14 - p^2 - 36 Cos[2 y] + p^2 Cos[2 y] + 12 p Sin[2 y]))]
         Sqrt[-11 + 12 p Tan[y] + 25 Tan[y]^2 - p^2 Tan[y]^2] + 
       Sqrt[2] p ArcTanh[(-11 - 
          6 I p + (-25 I + 6 p + I p^2) Tan[y])/((-6 I + 
            p) Sqrt[-11 + 12 p Tan[y] - (-25 + p^2) Tan[y]^2])] Cos[
         y] Sin[1] Sqrt[-((6 Cos[y] - p Sin[y])^2/(
         14 - p^2 - 36 Cos[2 y] + p^2 Cos[2 y] + 12 p Sin[2 y]))]
         Sqrt[-11 + 12 p Tan[y] + 25 Tan[y]^2 - p^2 Tan[y]^2] - 
       Sqrt[2] p ArcTanh[(
         11 - 6 I p + 
          I (-25 + 6 I p + p^2) Tan[y])/((6 I + p) Sqrt[-11 + 
           12 p Tan[y] - (-25 + p^2) Tan[y]^2])] Cos[y] Sin[
         1] Sqrt[-((6 Cos[y] - p Sin[y])^2/(
         14 - p^2 - 36 Cos[2 y] + p^2 Cos[2 y] + 12 p Sin[2 y]))]
         Sqrt[-11 + 12 p Tan[y] + 25 Tan[y]^2 - 
         p^2 Tan[y]^2])/(2 (6 Cos[1] - p Sin[1]) (-6 Cos[y] + 
         p Sin[y]))]}}*)

Visualization

cp1 = Plot[Evaluate[t[y] /. sol1 /. p -> 9.8], {y, -Pi, Pi}, 
  FrameLabel -> {"y", "t"}, PlotStyle -> Red, Frame -> True, 
  PlotLabel -> Row[{"E = ", E0}]]

fig2 = Show[
  Show[MapAt[Translate[#, {0, 2 Pi}] &, cp1, 1], cp1, 
   MapAt[Translate[#, {0, -2 Pi}] &, cp1, 1], PlotRange -> All], 
  Graphics[{Green, Thick, Line[{{0, -8}, {0, 13}}], 
    Line[{{Pi, -8}, {Pi, 13}}]}]]

Figure 2

Update 1 We can compare numerical solution proposed by Ulrich Neumann with exact solution as follows

erg[p_?NumericQ, E0_?NumericQ, M_?NumericQ] := 
 NDSolve[{t'[\[Sigma]]^2 == (E0*Cos[\[Sigma]] - 
        p*Sin[\[Sigma]])^2/((E0*Cos[\[Sigma]] - p*Sin[\[Sigma]])^2 - 
       M^2), t[1] == 2}, t, {\[Sigma], -Pi, Pi}, 
  Method -> {"StiffnessSwitching", 
    Method -> {"ExplicitRungeKutta", Automatic}}, 
  WorkingPrecision -> 30]

Visualization in a case of E0=6 with cp1 computed above

Show[cp1,Plot[Evaluate[Chop[t[y] /. erg[4911/500, 6, 5]]], {y, -Pi, Pi}, 
  PlotStyle -> {{Blue, Dashed}, {Blue, Dashed}}] // Quiet]

Figure 3

$\endgroup$
2
  • $\begingroup$ Nice solution too (+1)! Unfortunately Mathematica v12.2 can't reproduce sol1, that's why I tried numerical solution $\endgroup$ Commented Sep 16, 2023 at 7:09
  • $\begingroup$ @UlrichNeumann Oh, sorry, I forgot to add exact solution. $\endgroup$ Commented Sep 16, 2023 at 9:05
3
$\begingroup$

Try numerical solution NDSolve[...,Method->"StiffnessSwitching"].

Solutions seem to be complex, is this intended?

erg[p_?NumericQ, E0_?NumericQ, M_?NumericQ]  := 
NDSolve [{t'[\[Sigma]]^2 == (E0*Cos[\[Sigma]] -p*Sin[\[Sigma]])^2/((E0*Cos[\[Sigma]] - p*Sin[\[Sigma]])^2 - M^2), t[1] == 2}, t, {\[Sigma], -Pi, Pi}, 
Method -> "StiffnessSwitching"]

Plot[Evaluate[ReIm[t[y] /. erg[5.942, 0, 5]]], {y, -Pi, Pi}]

enter image description here

Plot[Evaluate[ReIm[t[y] /. erg[5.942, 6, 5]]], {y, -Pi, Pi}]

enter image description here

$\endgroup$
3
  • $\begingroup$ The solutions being complex could be the problem, because they should belong to the real domain. Is there a way to only search for real solutions? $\endgroup$
    – MarcelRomp
    Commented Sep 15, 2023 at 12:56
  • $\begingroup$ I don't know how to force real solutions. But have a look at your ode, t'[\[Sigma]]^2 sometimes is negativ. Check your ode! $\endgroup$ Commented Sep 15, 2023 at 13:06
  • $\begingroup$ @UlrichNeumann This is nice solution (+1). There are real solutions computed with special choice of E, M, p. See my answer. $\endgroup$ Commented Sep 16, 2023 at 4:18

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.