2
$\begingroup$ 
Clear["`.*"]
e = 0.92; sm = 15.; phi = 1.122; sii = ArcCos[e Cos[phi]];
NDSolve[
  {SI'[s] == Sin[PH[s]], PH'[s] == Sin[SI[s]]/e, 
   Y'[s] == Sin[PH[s]], X'[s] == Cos[PH[s]], PH[0] == phi, 
   SI[0] == sii, Y[0] == 0, X[0] == 0 }, 
  {SI, PH, Y, X}, {s, 0, sm}];
{si[u_], ph[u_], x[u_], y[u_]} = {SI[u], PH[u], X[u], Y[u]} /. First[%];
ParametricPlot[{x[s], y[s]}, {s, 0, sm}]

After a numerical solution is obtained as above, how can I convert it to an analytic solution?

EDIT1:

$$ \frac{\cos \psi}{\cos \phi } = e= \frac{\sin \psi\,\cdot \psi^{'}}{\sin \phi \,\cdot\phi^{'} }; $$

$$ \psi^{'} =\sin \phi, \quad \phi^{'} = \sin \psi /e ;$$

Conics_DEto Analyt_Soln

Conics tangent/directrix/eccentricity property above is used to set up first degree differential equations aimed to seek analytic solution from numerical (when axes parallel to coordinate axes) in the plane using Mathematica.

Improve this question
$\endgroup$
7
  • 1
    $\begingroup$ If I understand right, you cannot compute an exact solution from an approximate one. Or do you mean something else? $\endgroup$ – Michael E2 Jun 5 '17 at 21:13
  • 1
    $\begingroup$ Your plotted solution looks like an ellipse. Try fitting it to one. $\endgroup$ – bbgodfrey Jun 5 '17 at 21:47
  • 2
    $\begingroup$ 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:59
  • $\begingroup$ Also, Cos[SI[s]] == -C[1] + e Cos[PH[s]]. $\endgroup$ – bbgodfrey Jun 5 '17 at 23:51
  • 1
    $\begingroup$ Remember that the symbolic solutions involving elliptic functions that we are discussing provide s as a function of PH, not the inverse. Perhaps, it is possible also to obtain X and Y as solutions involving elliptic functions, but as functions of PH, not s. 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
4
$\begingroup$

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

enter image description here

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]

enter image description here

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

enter image description here

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]

enter image description here

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks soo much for the elaborate reply.. which serves me also as tutorial with application to conics, It should be also of general interest. Shall edit question reg conics differential equation background. $\endgroup$ – Narasimham Jun 12 '17 at 5:10
  • $\begingroup$ [enter link description here][2] [2]: math.stackexchange.com/questions/2319252/… // Gave a link to a relevant question, hope OK, $\endgroup$ – Narasimham Jun 12 '17 at 10:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.