Initial conditions with a non-linear ODE in Mathematica

Question

I'm trying to use Mathematica's NDSolve[] to compute a geodesic along a sphere using the coupled ODE:

x" - (x" . x) x = 0

The problem is that I can only enter initial conditions for x(0) and x'(0) and the solver is happy with the solution where x" = 0. The problem is that my geodesic on the sphere has the initial condition that x"(0) = -x(0), which I have no idea how to tell mathematica. If I add this as a condition, it says I'm adding True to the list of conditions.

Here is my code:

s1 = NDSolve[{x1''[t] - (x1[t] * x1''[t] + x2[t] * x2''[t] + x3[t]*x3''[t]) * x1[t] == 0, x2''[t] -  (x1[t] * x1''[t] + x2[t] * x2''[t] + x3[t]*x3''[t]) * x2[t] == 0, x3''[t] - (x1[t] * x1''[t] + x2[t] * x2''[t] + x3[t]*x3''[t]) * x3[t] == 0, x1[0] == 1, x2[0] == 0, x3[0] == 0, x1'[0] == 0, x2'[0] == 0, x3'[0] == 1} , { x1, x2, x3}, {t, -1, 1}][[1]]

I would like to modify this so that the initial acceleration is not zero but -x(0).

Thanks

Answer 1

I fixed this problem through a mathematical rearrangement rather than addressing my original issue:

Let V(t) be a vector field along x(t).

x . V = 0 implies d/dt (x . V) = (x' . V) + (x . V') = 0

So the equation D/dt V = V' - (x . V') x = V' + (x' . V) x holds This means the geodesic equation becomes: x" + (x' . x') x = 0 and so it can be solved using the initial conditions I originally had.

Thanks a lot Janus for going through and pointing out the various problems I was having including horrible code layout, I learnt a lot through your re-writing as well.

Answer 2

Well, as the error message says -- NDSolve only accepts initial conditions for derivatives of orders strictly less than the maximal order appearing in the ODE.
I have a feeling this is more of a mathematics question. Mathematically, {x''[0]=-x0, x[0]==x0}, doesn't define a unique solution - you'd have to do something along the lines of {x0.x''[0]==-1, x[0]==x0, x'[0]-x0 x0.x'[0]==v0} for that to work out (NDSolve would still fail with the same error). You do realize you will just get a great circle on the unit sphere, right?

By the way, here is how I would have coded up your example:

x[t_] = Table[Subscript[x, j][t], {j, 3}];
s1 = NDSolve[Flatten[Thread /@ #] &@{
       x''[t] - (x''[t].x[t]) x[t] == {0, 0, 0},
       x[0] == {1, 0, 0}, 
       x'[0] == {0, 0, 1}
     }, x[t], {t, -1, 1}]