Using spherical derivatives with NDSolve

Question

I am trying to plot the path of an object around a centre of mass. Obviously this needs spherical coordinates. The problem I encounter is with the derivatives. The Derivatives in spherical coordinates are for

(radial, azimuthal, polar) (r,theta,phi):

dr = r'
dϕ = r ϕ'
dθ = r Sin[ϕ] θ

Where the accelerations in theta and phi remain 0 for now.

But I need values for the initial velocity, this means that when I add the sinus to the NDSolve, it can't solve the set of equations; This is the furthest I've come. In truth, I'm not sure if this is the right approach or not, but solving the equations in Cartesian coordinates brings a lot more problems.

b = 20;
ball = First@NDSolve[{r''[t] == -9.81/r[t]^2, r[0] == 5, r'[0] == 0, 
  WhenEvent[r[t] == 2.25, r'[t] -> -0.9 r'[t]], θ''[t] == 0,
  θ[0] == 0,r[0] θ'[0] == π/4, ϕ''[t] == 0,
  ϕ[0] == π,r[0] ϕ'[0] == π/4}, {r, θ, ϕ}, {t, 0, b}];

The following doesn't work (third row, added sinus):

ball = First @ NDSolve[{r''[t] == -9.81/r[t]^2, r[0] == 5, r'[0] == 0, 
 WhenEvent[r[t] == 2.25, r'[t] -> -0.9 r'[t]], θ''[t] == 0,
 θ[0] == 0,r[0] Sin[ϕ[0]] θ'[0] == π/4, ϕ''[t] == 0,
 ϕ[0] == π,r[0] ϕ'[0] == π/4}, {r, θ, ϕ}, {t, 0, b}];

I then plot with:

Manipulate[
  Show[
    ParametricPlot3D[{r[t] Sin[θ[t]] Cos[ϕ[t]] /. ball, 
      r[t] Sin[θ[t]] Sin[ϕ[t]] /. ball, 
      r[t] Cos[θ[t]] /. ball}, {t, 0, a}, 
      PlotStyle -> {Gray, Dashed, Thick}, 
      PlotRange -> {{-5, 5}, {-5, 5}, {-5, 5}}], 
    ParametricPlot3D[{r[t] Sin[θ[t]] Cos[ϕ[t]] /. ball, 
      r[t] Sin[θ[t]] Sin[ϕ[t]] /. ball, 
      r[t] Cos[θ[t]] /. ball}, {t, a - 0.001, a}, 
      PlotStyle -> {Black, Thickness[0.05]}, 
      PlotRange -> {{-5, 5}, {-5, 5}, {-5, 5}}], 
    Graphics3D[Sphere[{0, 0, 0}, 2]]],
  {{a, b/100, "Time"}, b/100, b}]

Trying:

ball = First @ NDSolve[{x''[t]^2 + y''[t]^2 + z''[t]^2 ==
 -9.81^2/(x''[t]^2 + y''[t]^2 + z''[t]^2)^2, 
 x[0] == 3, y[0] == 3, z[0] == 3, x'[0] == 1, y'[0] == 0, 
 z'[0] == 0, 
 WhenEvent[x''[t]^2 + y''[t]^2 + z''[t]^2 == 2.25^2, 
  x''[t]^2 + y''[t]^2 + z''[t]^2 -> -0.9 x''[t]^2 + y''[t]^2 + 
  z''[t]^2]}, {x, y, z}, {t, 0, b}];

Returns an under-determined error.

I get that I have only one second derivative, but going all out (added row 3):

ball = First@NDSolve[{x''[t]^2 + y''[t]^2 + 
 z''[t]^2 == -9.81^2/(x''[t]^2 + y''[t]^2 + z''[t]^2)^2, 
 ArcTan[x[t], y[t]] == 0, ArcTan[z[t], Sqrt[x[t]^2 + y[t]^2]] == 0,
 x[0] == 3, y[0] == 3, z[0] == 3, x'[0] == 1, y'[0] == 0, 
 z'[0] == 0, 
 WhenEvent[x''[t]^2 + y''[t]^2 + z''[t]^2 == 2.25^2, 
  x''[t]^2 + y''[t]^2 + z''[t]^2 -> -0.9 x''[t]^2 + y''[t]^2 + 
  z''[t]^2]}, {x, y, z}, {t, 0, b}];

produces the messages:

Warning: the rule (x'')[t]^2+(y'')[t]^2+(z'')[t]^2->-0.9 (x'')[t]^2+(y'')[t]^2+(z''[t]^2 will not directly set the state because the left-hand side is not a list of state variables. >>

Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions. >>

Answer 1

Here is a systematic approach to this problem, starting from the Lagrangian: I denote the initial radial distance by r0 = 5. and use this value in writing the initial values for the other variables. The time of the simulation is tmax = 30.

Everything starts by defining the radius vector rVec in terms of spherical coordinates. Inserting this into the kinetic energy and subtracting the potential energy, we get the Lagrangian. The equations of motion eqs are automatically generated by using the function EulerEquations from the VariationalMethods package which I load at the beginning. That way, we can be sure that the translation to spherical coordinates is done correctly.

<< VariationalMethods`

rVec[t_] := 
 r[t] {Cos[ϕ[t]] Sin[θ[t]], 
   Sin[ϕ[t]] Sin[θ[t]], Cos[θ[t]]}

lagrangian = 1/2 Simplify[rVec'[t].rVec'[t]] + 9.81/r[t]

(*
==> 9.81/r[t] + 
 1/2 (Derivative[1][r][t]^2 + 
    r[t]^2 (Derivative[1][θ][t]^2 + 
       Sin[θ[t]]^2 Derivative[1][ϕ][t]^2))
*)

eqs = 
 EulerEquations[lagrangian, {r[t], θ[t], ϕ[t]}, t]

(*
==> {-(9.81/r[t]^2) + 
   r[t] (Derivative[1][θ][t]^2 + 
      Sin[θ[t]]^2 Derivative[1][ϕ][t]^2) - (
    r^′′)[t] == 0, 
 r[t] (-2 Derivative[1][r][t] Derivative[1][θ][t] + 
     r[t] (Cos[θ[t]] Sin[θ[t]] Derivative[1][ϕ][
          t]^2 - (θ^′′)[t])) == 
  0, -r[t] Sin[θ[
     t]] (2 Sin[θ[t]] Derivative[1][r][t] Derivative[
       1][ϕ][t] + 
     r[t] (2 Cos[θ[t]] Derivative[1][θ][t] Derivative[
          1][ϕ][t] + 
        Sin[θ[t]] (ϕ^′′)[t])) == 0}
*)

r0 = 5.;
tmax = 30;

sol = 
  First@NDSolve[
    Join[eqs, {r[0] == r0, 
      r'[0] == 0, θ[0] == 0.01, θ'[0] == Pi/4/r0, ϕ[0] == Pi, 
        ϕ'[0] == Pi/4/r0, 
      WhenEvent[r[t] == 2.25, r'[t] -> -0.9 r'[t]]}], {r[t], θ[
      t], ϕ[t]}, {t, 0, tmax}];

Show[Graphics3D[{Opacity[.5], Sphere[{0, 0, 0}, 2.25]}], 
 ParametricPlot3D[rVec[t] /. sol, {t, 0, tmax}, PlotPoints -> 100]]

In this plot, I had to increase the number of plot points to show the details of the bounces correctly. The motion is completely confined to a two-dimensional plane due to the conservation of angular momentum. The sphere is made slightly opaque to show the orbit better when you rotate the graphic.

I also had to change the starting condition for the polar angle away from $\theta=0$ because the spherical coordinate system is singular at the pole. So I chose θ[0] == 0.01 above.