Numerical Results
The equations in the question can be solved symbolically, but in terms of PH
, not s
. To begin, it is helpful to plot the numerical solution calculated in the question.
ParametricPlot[{x[s], y[s]}, {s, 0, sm}]

which is an ellipse. Plotting numerical results for the four dependent variables as functions of s
also is useful.
Plot[{si[s], x[s], y[s], ph[s]}, {s, 0, sm}, PlotRange -> All]

y[s]
does not look like x[s]
rescaled by e
and offset by a constant.
Symbolic Solution
Clear[e, phi]; sii = ArcCos[e Cos[phi]];
ds = DSolve[{SI'[s] == Sin[PH[s]] , PH'[s] == Sin[SI[s]]/e}, {SI[s], PH[s]}, s] // Last
(* {SI[s] -> ArcCos[-C[1] + e Cos[InverseFunction[-((4 I Cos[#1/2]^2 Sqrt[(
1 + C[1] - e Cos[#1])/((1 - e + C[1]) (1 + Cos[#1]))]
Sqrt[(1 - C[1] + e Cos[#1])/((1 + e - C[1]) (1 + Cos[#1]))] EllipticF[
I ArcSinh[Sqrt[(1 + e + C[1])/(1 - e + C[1])] Tan[#1/2]],
(1 - 2 e + e^2 - C[1]^2)/(1 + 2 e + e^2 - C[1]^2)])/(
Sqrt[(1 + e + C[1])/(2 - 2 e + 2 C[1])] Sqrt[
2 - e^2 - 2 C[1]^2 + 4 e C[1] Cos[#1] - e^2 Cos[2 #1]])) &][s/e + C[2]]]],
PH[s] -> InverseFunction[-((4 I Cos[#1/2]^2 Sqrt[(
1 + C[1] - e Cos[#1])/((1 - e + C[1]) (1 + Cos[#1]))] Sqrt[(
1 - C[1] + e Cos[#1])/((1 + e - C[1]) (1 + Cos[#1]))] EllipticF[
I ArcSinh[Sqrt[(1 + e + C[1])/(1 - e + C[1])] Tan[#1/2]],
(1 - 2 e + e^2 - C[1]^2)/(1 + 2 e + e^2 - C[1]^2)])/(
Sqrt[(1 + e + C[1])/(2 - 2 e + 2 C[1])] Sqrt[
2 - e^2 - 2 C[1]^2 + 4 e C[1] Cos[#1] - e^2 Cos[2 #1]])) &][s/e + C[2]]} *)
(DSolve
yields another symbolic solution which can be obtained from the one above by replacing s
by -s
and SI
by -SI
.) Now, determine the two constants of motion.
Flatten@Solve[Map[Cos, (First[ds] /. (Last[ds] /. Rule[z1_, z2_] -> Rule[z2, z1]))
/. Rule -> Equal], C[1]]
(* {C[1] -> e Cos[PH[s]] - Cos[SI[s]]} *)
Applying the initial conditions from the question then gives
% /. s -> 0 /. { PH[0] -> phi, SI[0] -> sii}
(* {C[1] -> 0} *)
and, because C[1]
vanishes, Cos[SI[s]] == e Cos[PH[s]]
, as also is evident from
ds0 = ds /. % /. Function[z_] :> Function @@ {z}
(* {SI[s] -> ArcCos[e Cos[
InverseFunction[-((4 I Cos[#1/2]^2 Sqrt[(1 - e Cos[#1])/((1 - e) (1 + Cos[#1]))]
Sqrt[(1 + e Cos[#1])/((1 + e) (1 + Cos[#1]))] EllipticF[I ArcSinh[Sqrt[(1 + e)/(1 - e)]
Tan[#1/2]], (1 - 2 e + e^2)/(1 + 2 e + e^2)])/(
Sqrt[(1 + e)/(2 - 2 e)] Sqrt[2 - e^2 - e^2 Cos[2 #1]])) &][s/e + C[2]]]],
PH[s] ->
InverseFunction[-((4 I Cos[#1/2]^2 Sqrt[(1 - e Cos[#1])/((1 - e) (1 + Cos[#1]))]
Sqrt[(1 + e Cos[#1])/((1 + e) (1 + Cos[#1]))] EllipticF[I ArcSinh[Sqrt[(1 + e)/(1 - e)]
Tan[#1/2]], (1 - 2 e + e^2)/(1 + 2 e + e^2)])/(
Sqrt[(1 + e)/(2 - 2 e)] Sqrt[2 - e^2 - e^2 Cos[2 #1]])) &][s/e + C[2]]]] *)
The symbolic expression and numerical value for C[2]
are obtained simply by
PH[s] /. ds0 /. s -> 0;
FullSimplify[%[[0, 1]][phi], Pi/2 > phi > -Pi/2 && 1 > e > 0]
% /. {e -> 0.92 , phi -> 1.122} // Chop
(* -((2 I EllipticF[I ArcSinh[Sqrt[-((1 + e)/(-1 + e))] Tan[phi/2]],
(-1 + e)^2/(1 + e)^2])/(1 + e)) *)
(* 1.91604 *)
To proceed, it convenient to determine X
and Y
in terms of PH
, which is accomplished by transforming the independent variable in their ODEs to PH
and employing Integrate
to obtain solutions up to a constant.
xs = FullSimplify[Integrate[e Cos[PH]/Sqrt[1 - e^2 Cos[PH]^2], PH], 0 < e < 1]
(* ArcTanh[(Sqrt[2] e Sin[PH])/Sqrt[2 - e^2 - e^2 Cos[2 PH]]] *)
ys = Integrate[e Sin[PH]/Sqrt[1 - e^2 Cos[PH]^2], PH]
(* -ArcSin[e Cos[PH]] *)
These two expressions, plus Cos[SI[s]] == e Cos[PH[s]]
, obtained earlier, give
Plot[{ArcCos[e Cos[PH] /. e -> 0.92], Evaluate[(xs - (xs /. PH -> ph[0])) /. e -> 0.92],
Evaluate[(ys - (ys /. PH -> ph[0])) /. e -> 0.92]}, {PH, ph[0], ph[sm]}]

This plot is identical to what can be obtained for the numerical solutions by
ParametricPlot[{{ph[s], si[s]}, {ph[s], x[s]}, {ph[s], y[s]}}, {s, 0, sm},
PlotRange -> All, AspectRatio -> 1/GoldenRatio]
Thus, SI
, X
, and Y
have been obtained as symbolic functions of PH
, and PH
as a seemingly symbolic function of s
in ds0
. Unfortunately, InverseFunction
is not really symbolic, and moreover does not return numerical values for the full range of s
. So, the final step is to obtain s
as a function of PH
, which can be done by extracting the function inside InverseFunction
and applying it to PH
, yielding for s
.
e ((PH[s] /. ds0)[[0, 1]][PH] - C[2])
(* e (-C[2] - (4 I Cos[PH/2]^2 Sqrt[(1 - e Cos[PH])/((1 - e) (1 + Cos[PH]))]
Sqrt[(1 + e Cos[PH])/((1 + e) (1 + Cos[PH]))]
EllipticF[I ArcSinh[Sqrt[(1 + e)/(1 - e)] Tan[PH/2]], (1 - 2 e + e^2)/
(1 + 2 e + e^2)])/(Sqrt[(1 + e)/(2 - 2 e)] Sqrt[2 - e^2 - e^2 Cos[2 PH]])) *)
However, this expression is multivalued, so its period must be determined,
fun = % /. e -> .92 /. C[2] -> 1.9160429061618474`;
per = Re@N[fun /. PH -> Pi - 10^-7] - Re@N[fun /. PH -> Pi + 10^-7]
(* 8.7583 *)
Finally, this gives s
vs PH
.
ParametricPlot[{Piecewise[{{fun, PH < Pi}, {fun + per, PH < 3 Pi},
{fun + 2 per, PH < 5 Pi}}], PH}, {PH, ph[0], ph[15]},
AspectRatio -> 1/GoldenRatio, PlotStyle -> Red]

All four dependent functions now have been obtained parametrically in terms of s
. Although I undertook this solution as a challenge, realistically the numerical solution is more useful for most purposes.
DSolve[{SI'[s] == Sin[PH[s]] , PH'[s] == Sin[SI[s]]/e}, {SI[s], PH[s]}, s]
yields a symbolic solution, although not a very convenient one. $\endgroup$ – bbgodfrey Jun 5 '17 at 21:59Cos[SI[s]] == -C[1] + e Cos[PH[s]]
. $\endgroup$ – bbgodfrey Jun 5 '17 at 23:51s
as a function ofPH
, not the inverse. Perhaps, it is possible also to obtainX
andY
as solutions involving elliptic functions, but as functions ofPH
, nots
. The fundamental issue from my perspective is, what precisely are you seeking. To say that you want "an analytic solution" is not sufficient, unless you define precisely what you mean by "analytical". Elliptic functions, for instance, are analytical, but are integrals of them, or their inverses? $\endgroup$ – bbgodfrey Jun 6 '17 at 13:15