I would be extremely grateful for any help regarding the following code I wrote and the errors it produces. In this code I am investigating the behaviour of a massive particle trapped in the vicinity of the gravitational Bessel wave vertex.
I am really sorry for the lengthy expressions - I have no idea how to shorten it in any way that would help significantly. The table defined below is the right hand side of the geodesic equation. I have checked its correctness with my advisor, obtaining the same result as he did. The parameter "a" refers to the amplitude of the wave, while "b" and "c" are responsible for the wavevector quantities.
Remove["Global`*"]
coords = {t, x, y, z}
n = Length[coords];
geodesic = {BesselJ[3,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] (-((a b^4 Sqrt[b^2 + c^2] x (x^2 - 3 y^2) (x')^2)/(x^2 +
y^2)^(3/2)) + (2 a b^4 Sqrt[
b^2 + c^2] y (-3 x^2 + y^2) x' y')/(x^2 + y^2)^(3/
2) + (a b^4 Sqrt[b^2 + c^2] x (x^2 - 3 y^2) (y')^2)/(x^2 +
y^2)^(3/2)) +
Cos[Sqrt[b^2 + c^2] t -
c z] ((a b^4 Sqrt[b^2 + c^2] y (-3 x^2 + y^2) (x')^2)/(x^2 +
y^2)^(3/2) + (2 a b^4 Sqrt[
b^2 + c^2] x (x^2 - 3 y^2) x' y')/(x^2 + y^2)^(3/
2) - (a b^4 Sqrt[
b^2 + c^2] (-3 x^2 y + y^3) (y')^2)/(x^2 + y^2)^(3/2))) +
BesselJ[0,
b Sqrt[x^2 + y^2]] (4 a b Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^3 Cos[
Sqrt[b^2 + c^2] t - c z] x' z' +
4 a b Sqrt[b^2 + c^2] (c + Sqrt[b^2 + c^2])^3 Sin[
Sqrt[b^2 + c^2] t - c z] y' z') +
BesselJ[2,
b Sqrt[x^2 + y^2]] (Cos[
Sqrt[b^2 + c^2] t -
c z] ((4 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) (x^2 - y^2) x' z')/(x^2 +
y^2) + (8 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) x y y' z')/(x^2 + y^2)) -
Sin[
Sqrt[b^2 + c^2] t -
c z] (-((8 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) x y x' z')/(x^2 +
y^2)) - (4 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) (-x^2 +
y^2) y' z')/(x^2 + y^2))) +
BesselJ[1,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] ((a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x (x')^2)/
Sqrt[x^2 +
y^2] + (2 a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^4 y x' y')/
Sqrt[x^2 +
y^2] - (a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x (y')^2)/
Sqrt[x^2 +
y^2] + (4 a b^2 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 x (z')^2)/
Sqrt[x^2 + y^2]) +
Cos[Sqrt[b^2 + c^2] t -
c z] (-((a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) y (x')^2)/
Sqrt[x^2 + y^2]) + (2 a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^4 x x' y')/
Sqrt[x^2 +
y^2] + (a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) y (y')^2)/
Sqrt[
x^2 + y^2] + (4 a b^2 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 y (z')^2)/
Sqrt[x^2 + y^2])),
BesselJ[4,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] ((a b^5 x y (-3 x^2 +
y^2) (x')^2)/(2 (x^2 + y^2)^2) + (a b^5 y^2 (-3 x^2 +
y^2) x' y')/(x^2 +
y^2)^2 - (a b^5 x y (x^2 +
5 y^2) (y')^2)/(2 (x^2 + y^2)^2)) +
Cos[Sqrt[b^2 + c^2] t -
c z] ((a b^5 x^2 (x^2 -
3 y^2) (x')^2)/(2 (x^2 + y^2)^2) + (a b^5 x y (x^2 -
3 y^2) x' y')/(x^2 +
y^2)^2 + (a b^5 (x^4 + 3 x^2 y^2 -
2 y^4) (y')^2)/(2 (x^2 + y^2)^2))) +
BesselJ[3,
b Sqrt[x^2 + y^2]] (Cos[
Sqrt[b^2 + c^2] t -
c z] ((2 a b^4 Sqrt[b^2 + c^2] y (-3 x^2 + y^2) t' x')/(x^2 +
y^2)^(3/2) + (3 a b^4 y (-3 x^2 y + y^3) (x')^2)/(x^2 +
y^2)^(5/2) + (2 a b^4 Sqrt[
b^2 + c^2] x (x^2 - 3 y^2) t' y')/(x^2 + y^2)^(3/
2) - (6 a b^4 x y (-3 x^2 + y^2) x' y')/(x^2 + y^2)^(5/
2) + (3 a b^4 (-2 x^4 + 3 x^2 y^2 + y^4) (y')^2)/(x^2 +
y^2)^(5/2) - (2 a b^4 c y (-3 x^2 + y^2) x' z')/(x^2 +
y^2)^(3/
2) + (2 a b^4 x (4 c y^2 +
Sqrt[b^2 + c^2] (x^2 + y^2)) y' z')/(x^2 + y^2)^(3/2)) -
Sin[Sqrt[b^2 + c^2] t -
c z] (-((2 a b^4 Sqrt[b^2 + c^2] x (x^2 - 3 y^2) t' x')/(x^2 +
y^2)^(3/2)) - (3 a b^4 x y (x^2 - 3 y^2) (x')^2)/(x^2 +
y^2)^(5/2) + (2 a b^4 Sqrt[
b^2 + c^2] y (-3 x^2 + y^2) t' y')/(x^2 + y^2)^(3/
2) + (6 a b^4 x^2 (x^2 - 3 y^2) x' y')/(x^2 + y^2)^(5/
2) + (3 a b^4 x y (5 x^2 + y^2) (y')^2)/(x^2 + y^2)^(5/
2) + (2 a b^4 c x (x^2 - 3 y^2) x' z')/(x^2 + y^2)^(3/
2) - (2 a b^4 y (-2 c (x^2 - y^2) +
Sqrt[b^2 + c^2] (x^2 + y^2)) y' z')/(x^2 + y^2)^(3/2))) +
BesselJ[1,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] ((2 a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x t' x')/
Sqrt[x^2 +
y^2] - (a (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y (x')^2)/(x^2 + y^2)^(3/
2) + (2 a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^4 y t' y')/
Sqrt[
x^2 + y^2] + (2 a (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x^2 x' y')/(x^2 + y^2)^(3/
2) + (a (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y (y')^2)/(x^2 + y^2)^(3/
2) - (2 a c (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x x' z')/
Sqrt[x^2 +
y^2] - (2 a (c + Sqrt[b^2 + c^2]) (c^4 +
3 c^3 Sqrt[b^2 + c^2] - b^2 (b^2 + 2 (b^2 + c^2)) +
c^2 (-2 b^2 + 3 (b^2 + c^2)) +
c (-4 b^2 Sqrt[b^2 + c^2] + (b^2 + c^2)^(3/2))) y y' z')/
Sqrt[x^2 +
y^2] - (4 a b^2 (c + Sqrt[b^2 + c^2])^2 x y (z')^2)/(x^2 +
y^2)^(3/2)) +
Cos[Sqrt[b^2 + c^2] t -
c z] (-((2 a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) y t' x')/
Sqrt[x^2 + y^2]) + (a (c + Sqrt[b^2 + c^2])^2 (-b^2 +
2 c^2 + 2 c Sqrt[b^2 + c^2]) y^2 (x')^2)/(x^2 + y^2)^(3/
2) + (2 a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^4 x t' y')/
Sqrt[x^2 +
y^2] - (2 a (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y x' y')/(x^2 + y^2)^(3/
2) + (a (c +
Sqrt[b^2 + c^2])^2 (2 (b^2 + c^2) x^2 + (3 b^2 +
c^2) y^2 + c^2 (2 x^2 + y^2) +
2 c Sqrt[b^2 + c^2] (2 x^2 + y^2)) (y')^2)/(x^2 +
y^2)^(3/
2) + (2 a c (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) y x' z')/
Sqrt[x^2 +
y^2] - (2 a (c + Sqrt[b^2 + c^2]) (b^4 + c^4 +
3 c^3 Sqrt[b^2 + c^2] + 2 b^2 (b^2 + c^2) +
c^2 (2 b^2 + 3 (b^2 + c^2)) +
c (4 b^2 Sqrt[b^2 + c^2] + (b^2 + c^2)^(3/2))) x y' z')/
Sqrt[x^2 +
y^2] - (4 a b^2 (c + Sqrt[b^2 + c^2])^2 y^2 (z')^2)/(x^2 +
y^2)^(3/2))) +
BesselJ[2,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] (-((a b x y ((-3 b^4 +
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) x^2 + (2 b^2 +
c^2)^2 y^2 + c^4 (x^2 + y^2) +
4 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
4 c Sqrt[b^2 + c^2] (2 b^2 + c^2) (x^2 + y^2) +
2 c^2 (b^2 + 3 (b^2 + c^2)) (x^2 +
y^2)) (x')^2)/(2 (x^2 +
y^2)^2)) - (a b y^2 ((-3 b^4 +
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) x^2 + (2 b^2 +
c^2)^2 y^2 + c^4 (x^2 + y^2) +
4 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
4 c Sqrt[b^2 + c^2] (2 b^2 + c^2) (x^2 + y^2) +
2 c^2 (b^2 + 3 (b^2 + c^2)) (x^2 + y^2)) x' y')/(x^2 +
y^2)^2 + (a b x y ((2 b^2 + c^2)^2 x^2 + (5 b^4 +
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) y^2 +
c^4 (x^2 + y^2) + 4 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
4 c Sqrt[b^2 + c^2] (2 b^2 + c^2) (x^2 + y^2) +
2 c^2 (b^2 + 3 (b^2 + c^2)) (x^2 +
y^2)) (y')^2)/(2 (x^2 + y^2)^2) - (8 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) x y t' z')/(x^2 +
y^2) + (8 a b^3 (c + Sqrt[b^2 + c^2]) y y' z')/(x^2 +
y^2) - (2 a b^3 (-3 c + Sqrt[b^2 + c^2]) (c +
Sqrt[b^2 + c^2]) x y (z')^2)/(x^2 + y^2)) +
Cos[Sqrt[b^2 + c^2] t -
c z] (-((a b x^2 (c^4 x^2 + (-3 b^4 -
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) y^2 +
c^4 (x^2 + y^2) + 8 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) -
2 c^2 (b^2 - 3 (b^2 + c^2)) (x^2 +
y^2)) (x')^2)/(2 (x^2 +
y^2)^2)) - (a b x y (c^4 x^2 + (-3 b^4 -
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) y^2 +
c^4 (x^2 + y^2) + 8 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) -
2 c^2 (b^2 - 3 (b^2 + c^2)) (x^2 + y^2)) x' y')/(x^2 +
y^2)^2 -
1/(2 (x^2 + y^2)^2) a b ((2 b^2 + c^2)^2 x^4 + (3 b^4 +
2 b^2 (b^2 + c^2) + 3 (b^2 + c^2)^2) x^2 y^2 +
2 (-b^4 + (b^2 + c^2)^2) y^4 +
4 c Sqrt[
b^2 + c^2] (x^2 + y^2) ((2 b^2 + c^2) x^2 +
2 (b^2 + c^2) y^2) +
2 c^2 (x^2 + y^2) ((b^2 + 3 (b^2 + c^2)) x^2 +
6 (b^2 + c^2) y^2) + c^4 (x^4 + 3 x^2 y^2 + 2 y^4) +
4 c^3 Sqrt[
b^2 + c^2] (x^4 + 3 x^2 y^2 +
2 y^4)) (y')^2 + (4 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) (x^2 - y^2) t' z')/(x^2 +
y^2) - (8 a b^3 (c + Sqrt[b^2 + c^2]) x y' z')/(x^2 +
y^2) - (2 a b^3 (c +
Sqrt[b^2 + c^2]) (-Sqrt[b^2 + c^2] x^2 +
c (x^2 - 2 y^2)) (z')^2)/(x^2 + y^2))) +
BesselJ[0,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] ((a b (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y (x')^2)/(2 (x^2 +
y^2)) + (a b (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) y^2 x' y')/(x^2 +
y^2) - (a b (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y (y')^2)/(2 (x^2 +
y^2)) + (2 a b^3 (c +
Sqrt[b^2 + c^2])^2 x y (z')^2)/(x^2 + y^2)) +
Cos[Sqrt[b^2 + c^2] t -
c z] ((a b (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x^2 (x')^2)/(2 (x^2 +
y^2)) + (a b (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y x' y')/(x^2 +
y^2) + (a b (c + Sqrt[b^2 + c^2])^2 ((3 b^2 + c^2) x^2 +
2 (b^2 + c^2) y^2 + c^2 (x^2 + 2 y^2) +
2 c Sqrt[b^2 + c^2] (x^2 + 2 y^2)) (y')^2)/(2 (x^2 +
y^2)) +
4 a b Sqrt[
b^2 + c^2] (c +
Sqrt[b^2 + c^2])^3 t' z' - (2 a b (c +
Sqrt[b^2 + c^2])^2 (b^2 x^2 + 2 c^2 (x^2 + y^2) +
2 c Sqrt[b^2 + c^2] (x^2 + y^2)) (z')^2)/(x^2 + y^2))),
BesselJ[4,
b Sqrt[x^2 + y^2]] (Cos[
Sqrt[b^2 + c^2] t -
c z] ((a b^5 x y (5 x^2 +
y^2) (x')^2)/(2 (x^2 + y^2)^2) - (a b^5 x^2 (x^2 -
3 y^2) x' y')/(x^2 +
y^2)^2 - (a b^5 x y (x^2 -
3 y^2) (y')^2)/(2 (x^2 + y^2)^2)) -
Sin[Sqrt[b^2 + c^2] t -
c z] (-((a b^5 (-2 x^4 + 3 x^2 y^2 +
y^4) (x')^2)/(2 (x^2 + y^2)^2)) - (a b^5 x y (-3 x^2 +
y^2) x' y')/(x^2 +
y^2)^2 - (a b^5 y^2 (-3 x^2 +
y^2) (y')^2)/(2 (x^2 + y^2)^2))) +
BesselJ[3,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] ((2 a b^4 Sqrt[b^2 + c^2] y (-3 x^2 + y^2) t' x')/(x^2 +
y^2)^(3/
2) - (3 a b^4 (x^4 + 3 x^2 y^2 - 2 y^4) (x')^2)/(x^2 +
y^2)^(5/2) + (2 a b^4 Sqrt[
b^2 + c^2] x (x^2 - 3 y^2) t' y')/(x^2 + y^2)^(3/
2) + (6 a b^4 x y (x^2 - 3 y^2) x' y')/(x^2 + y^2)^(5/
2) - (3 a b^4 x^2 (x^2 - 3 y^2) (y')^2)/(x^2 + y^2)^(5/
2) + (2 a b^4 y (4 c x^2 +
Sqrt[b^2 + c^2] (x^2 + y^2)) x' z')/(x^2 + y^2)^(3/
2) - (2 a b^4 c x (x^2 - 3 y^2) y' z')/(x^2 + y^2)^(3/
2)) + Cos[
Sqrt[b^2 + c^2] t -
c z] ((2 a b^4 Sqrt[b^2 + c^2] x (x^2 - 3 y^2) t' x')/(x^2 +
y^2)^(3/2) - (3 a b^4 x y (x^2 + 5 y^2) (x')^2)/(x^2 +
y^2)^(5/2) - (2 a b^4 Sqrt[
b^2 + c^2] (-3 x^2 y + y^3) t' y')/(x^2 + y^2)^(3/
2) - (6 a b^4 y (-3 x^2 y + y^3) x' y')/(x^2 + y^2)^(5/
2) + (3 a b^4 x y (-3 x^2 + y^2) (y')^2)/(x^2 + y^2)^(5/
2) - (2 a b^4 x (2 c (x^2 - y^2) +
Sqrt[b^2 + c^2] (x^2 + y^2)) x' z')/(x^2 + y^2)^(3/
2) + (2 a b^4 c (-3 x^2 y + y^3) y' z')/(x^2 + y^2)^(3/
2))) + BesselJ[1,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] ((2 a Sqrt[b^2 + c^2] (c + Sqrt[b^2 + c^2])^4 y t' x')/
Sqrt[x^2 +
y^2] - (a (c + Sqrt[b^2 + c^2])^2 ((3 b^2 + c^2) x^2 +
2 (b^2 + c^2) y^2 + c^2 (x^2 + 2 y^2) +
2 c Sqrt[b^2 + c^2] (x^2 + 2 y^2)) (x')^2)/(x^2 +
y^2)^(3/2) - (2 a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x t' y')/
Sqrt[x^2 +
y^2] + (2 a (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y x' y')/(x^2 + y^2)^(3/
2) - (a (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x^2 (y')^2)/(x^2 + y^2)^(3/
2) - (2 a (c + Sqrt[b^2 + c^2]) (b^4 + c^4 +
3 c^3 Sqrt[b^2 + c^2] + 2 b^2 (b^2 + c^2) +
c^2 (2 b^2 + 3 (b^2 + c^2)) +
c (4 b^2 Sqrt[b^2 + c^2] + (b^2 + c^2)^(3/2))) y x' z')/
Sqrt[x^2 +
y^2] + (2 a c (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y' z')/
Sqrt[x^2 +
y^2] + (4 a b^2 (c + Sqrt[b^2 + c^2])^2 x^2 (z')^2)/(x^2 +
y^2)^(3/2)) +
Cos[Sqrt[b^2 + c^2] t -
c z] ((2 a Sqrt[b^2 + c^2] (c + Sqrt[b^2 + c^2])^4 x t' x')/
Sqrt[x^2 +
y^2] - (a (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y (x')^2)/(x^2 + y^2)^(3/
2) + (2 a Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) y t' y')/
Sqrt[x^2 +
y^2] - (2 a (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) y^2 x' y')/(x^2 + y^2)^(3/
2) + (a (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y (y')^2)/(x^2 + y^2)^(3/
2) - (2 a (c + Sqrt[b^2 + c^2]) (c^4 +
3 c^3 Sqrt[b^2 + c^2] - b^2 (b^2 + 2 (b^2 + c^2)) +
c^2 (-2 b^2 + 3 (b^2 + c^2)) +
c (-4 b^2 Sqrt[b^2 + c^2] + (b^2 + c^2)^(3/2))) x x' z')/
Sqrt[x^2 +
y^2] - (2 a c (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) y y' z')/
Sqrt[x^2 +
y^2] + (4 a b^2 (c + Sqrt[b^2 + c^2])^2 x y (z')^2)/(x^2 +
y^2)^(3/2))) +
BesselJ[2,
b Sqrt[x^2 + y^2]] (Cos[
Sqrt[b^2 + c^2] t -
c z] (-((a b x y ((5 b^4 +
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) x^2 + (2 b^2 +
c^2)^2 y^2 + c^4 (x^2 + y^2) +
4 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
4 c Sqrt[b^2 + c^2] (2 b^2 + c^2) (x^2 + y^2) +
2 c^2 (b^2 + 3 (b^2 + c^2)) (x^2 +
y^2)) (x')^2)/(2 (x^2 +
y^2)^2)) + (a b x^2 ((2 b^2 + c^2)^2 x^2 + (-3 b^4 +
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) y^2 +
c^4 (x^2 + y^2) + 4 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
4 c Sqrt[b^2 + c^2] (2 b^2 + c^2) (x^2 + y^2) +
2 c^2 (b^2 + 3 (b^2 + c^2)) (x^2 + y^2)) x' y')/(x^2 +
y^2)^2 + (a b x y ((2 b^2 + c^2)^2 x^2 + (-3 b^4 +
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) y^2 +
c^4 (x^2 + y^2) + 4 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
4 c Sqrt[b^2 + c^2] (2 b^2 + c^2) (x^2 + y^2) +
2 c^2 (b^2 + 3 (b^2 + c^2)) (x^2 +
y^2)) (y')^2)/(2 (x^2 + y^2)^2) + (8 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) x y t' z')/(x^2 +
y^2) + (8 a b^3 (c + Sqrt[b^2 + c^2]) x x' z')/(x^2 +
y^2) + (2 a b^3 (-3 c + Sqrt[b^2 + c^2]) (c +
Sqrt[b^2 + c^2]) x y (z')^2)/(x^2 + y^2)) -
Sin[Sqrt[b^2 + c^2] t -
c z] (1/(2 (x^2 +
y^2)^2) a b (-2 (b^4 - (b^2 + c^2)^2) x^4 + (3 b^4 +
2 b^2 (b^2 + c^2) +
3 (b^2 + c^2)^2) x^2 y^2 + (2 b^2 + c^2)^2 y^4 +
4 c Sqrt[
b^2 + c^2] (x^2 +
y^2) (2 (b^2 + c^2) x^2 + (2 b^2 + c^2) y^2) +
2 c^2 (x^2 +
y^2) (6 (b^2 + c^2) x^2 + (b^2 + 3 (b^2 + c^2)) y^2) +
c^4 (2 x^4 + 3 x^2 y^2 + y^4) +
4 c^3 Sqrt[
b^2 + c^2] (2 x^4 + 3 x^2 y^2 +
y^4)) (x')^2 + (a b x y ((-3 b^4 -
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) x^2 + c^4 y^2 +
c^4 (x^2 + y^2) + 8 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) -
2 c^2 (b^2 - 3 (b^2 + c^2)) (x^2 + y^2)) x' y')/(x^2 +
y^2)^2 + (a b y^2 ((-3 b^4 -
2 b^2 (b^2 + c^2) + (b^2 + c^2)^2) x^2 + c^4 y^2 +
c^4 (x^2 + y^2) + 8 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) -
2 c^2 (b^2 - 3 (b^2 + c^2)) (x^2 +
y^2)) (y')^2)/(2 (x^2 + y^2)^2) - (4 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) (-x^2 +
y^2) t' z')/(x^2 +
y^2) - (8 a b^3 (c + Sqrt[b^2 + c^2]) y x' z')/(x^2 +
y^2) + (2 a b^3 (c +
Sqrt[b^2 + c^2]) (-Sqrt[b^2 + c^2] y^2 +
c (-2 x^2 + y^2)) (z')^2)/(x^2 + y^2))) +
BesselJ[0,
b Sqrt[x^2 + y^2]] (Cos[
Sqrt[b^2 + c^2] t -
c z] ((a b (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y (x')^2)/(2 (x^2 +
y^2)) - (a b (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x^2 x' y')/(x^2 +
y^2) - (a b (c + Sqrt[b^2 + c^2])^2 (3 b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y (y')^2)/(2 (x^2 +
y^2)) - (2 a b^3 (c +
Sqrt[b^2 + c^2])^2 x y (z')^2)/(x^2 + y^2)) -
Sin[Sqrt[b^2 + c^2] t -
c z] (-((a b (c +
Sqrt[b^2 + c^2])^2 (2 (b^2 + c^2) x^2 + (3 b^2 +
c^2) y^2 + c^2 (2 x^2 + y^2) +
2 c Sqrt[b^2 + c^2] (2 x^2 + y^2)) (x')^2)/(2 (x^2 +
y^2))) - (a b (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) x y x' y')/(x^2 +
y^2) - (a b (c + Sqrt[b^2 + c^2])^2 (-b^2 + 2 c^2 +
2 c Sqrt[b^2 + c^2]) y^2 (y')^2)/(2 (x^2 + y^2)) -
4 a b Sqrt[
b^2 + c^2] (c +
Sqrt[b^2 + c^2])^3 t' z' + (2 a b (c +
Sqrt[b^2 + c^2])^2 (b^2 y^2 + 2 c^2 (x^2 + y^2) +
2 c Sqrt[b^2 + c^2] (x^2 + y^2)) (z')^2)/(x^2 + y^2))),
BesselJ[3,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] ((a b^4 x (2 Sqrt[b^2 + c^2] (x^2 - y^2) +
c (x^2 + y^2)) (x')^2)/(x^2 + y^2)^(3/2) - (2 a b^4 Sqrt[
b^2 + c^2] y (-3 x^2 + y^2) x' y')/(x^2 + y^2)^(3/
2) + (a b^4 x (4 Sqrt[b^2 + c^2] y^2 +
c (x^2 + y^2)) (y')^2)/(x^2 + y^2)^(3/2)) +
Cos[Sqrt[b^2 + c^2] t -
c z] ((a b^4 y (4 Sqrt[b^2 + c^2] x^2 +
c (x^2 + y^2)) (x')^2)/(x^2 + y^2)^(3/2) - (2 a b^4 Sqrt[
b^2 + c^2] x (x^2 - 3 y^2) x' y')/(x^2 + y^2)^(3/
2) + (a b^4 y (2 Sqrt[b^2 + c^2] (-x^2 + y^2) +
c (x^2 + y^2)) (y')^2)/(x^2 + y^2)^(3/2))) +
BesselJ[0,
b Sqrt[x^2 + y^2]] (Cos[
Sqrt[b^2 + c^2] t -
c z] (4 a b Sqrt[
b^2 + c^2] (c +
Sqrt[b^2 + c^2])^3 t' x' + (4 a b^3 (c +
Sqrt[b^2 + c^2])^2 x^2 x' z')/(x^2 +
y^2) + (4 a b^3 (c + Sqrt[b^2 + c^2])^2 x y y' z')/(x^2 +
y^2)) -
Sin[Sqrt[b^2 + c^2] t -
c z] (-4 a b Sqrt[
b^2 + c^2] (c +
Sqrt[b^2 + c^2])^3 t' y' - (4 a b^3 (c +
Sqrt[b^2 + c^2])^2 x y x' z')/(x^2 +
y^2) - (4 a b^3 (c + Sqrt[b^2 + c^2])^2 y^2 y' z')/(x^2 +
y^2))) +
BesselJ[2,
b Sqrt[x^2 + y^2]] (Cos[
Sqrt[b^2 + c^2] t -
c z] ((4 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) (x^2 - y^2) t' x')/(x^2 +
y^2) - (8 a b^3 (c + Sqrt[b^2 + c^2]) y (-x^2 +
y^2) (x')^2)/(x^2 + y^2)^2 + (8 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) x y t' y')/(x^2 +
y^2) - (8 a b^3 (c + Sqrt[b^2 + c^2]) x (x^2 -
3 y^2) x' y')/(x^2 +
y^2)^2 - (16 a b^3 (c +
Sqrt[b^2 + c^2]) x^2 y (y')^2)/(x^2 +
y^2)^2 - (4 a b^3 (c +
Sqrt[b^2 + c^2])^2 x^2 x' z')/(x^2 +
y^2) - (4 a b^3 (c + Sqrt[b^2 + c^2])^2 x y y' z')/(x^2 +
y^2)) -
Sin[Sqrt[b^2 + c^2] t -
c z] (-((8 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) x y t' x')/(x^2 +
y^2)) - (16 a b^3 (c +
Sqrt[b^2 + c^2]) x y^2 (x')^2)/(x^2 +
y^2)^2 - (4 a b^3 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2]) (-x^2 +
y^2) t' y')/(x^2 +
y^2) - (8 a b^3 (c + Sqrt[b^2 + c^2]) y (-3 x^2 +
y^2) x' y')/(x^2 +
y^2)^2 + (8 a b^3 (c + Sqrt[b^2 + c^2]) x (-x^2 +
y^2) (y')^2)/(x^2 +
y^2)^2 + (4 a b^3 (c +
Sqrt[b^2 + c^2])^2 x y x' z')/(x^2 +
y^2) + (4 a b^3 (c + Sqrt[b^2 + c^2])^2 y^2 y' z')/(x^2 +
y^2))) +
BesselJ[1,
b Sqrt[x^2 +
y^2]] (-Sin[
Sqrt[b^2 + c^2] t -
c z] ((a (c + Sqrt[b^2 + c^2]) x (c^4 (x^2 + y^2) +
3 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
c Sqrt[b^2 + c^2] (7 b^2 + c^2) (x^2 + y^2) +
c^2 (2 b^2 + 3 (b^2 + c^2)) (x^2 + y^2) +
2 b^2 (-(b^2 - 2 (b^2 + c^2)) x^2 + (b^2 +
2 (b^2 + c^2)) y^2)) (x')^2)/(x^2 + y^2)^(3/
2) + (2 a (c + Sqrt[b^2 + c^2]) y (c^4 (x^2 + y^2) +
3 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
c Sqrt[b^2 + c^2] (5 b^2 + c^2) (x^2 + y^2) +
c^2 (2 b^2 + 3 (b^2 + c^2)) (x^2 + y^2) +
b^2 ((-3 b^2 + 2 (b^2 + c^2)) x^2 + (b^2 +
2 (b^2 + c^2)) y^2)) x' y')/(x^2 + y^2)^(3/
2) - (a (c + Sqrt[b^2 + c^2]) x (4 b^4 y^2 +
c^4 (x^2 + y^2) + 3 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
c Sqrt[b^2 + c^2] (3 b^2 + c^2) (x^2 + y^2) +
c^2 (2 b^2 + 3 (b^2 + c^2)) (x^2 + y^2)) (y')^2)/(x^2 +
y^2)^(3/2) + (8 a b^2 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 x t' z')/
Sqrt[x^2 +
y^2] + (8 a b^2 (c + Sqrt[b^2 + c^2])^2 x y x' z')/(x^2 +
y^2)^(3/
2) - (8 a b^2 (c + Sqrt[b^2 + c^2])^2 x^2 y' z')/(x^2 +
y^2)^(3/2) - (4 a b^2 c (c + Sqrt[b^2 + c^2])^2 x (z')^2)/
Sqrt[x^2 + y^2]) +
Cos[Sqrt[b^2 + c^2] t -
c z] (-((a (c + Sqrt[b^2 + c^2]) y (4 b^4 x^2 +
c^4 (x^2 + y^2) + 3 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
c Sqrt[b^2 + c^2] (3 b^2 + c^2) (x^2 + y^2) +
c^2 (2 b^2 + 3 (b^2 + c^2)) (x^2 + y^2)) (x')^2)/(x^2 +
y^2)^(3/2)) + (2 a (c +
Sqrt[b^2 + c^2]) x (c^4 (x^2 + y^2) +
3 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
c Sqrt[b^2 + c^2] (5 b^2 + c^2) (x^2 + y^2) +
c^2 (2 b^2 + 3 (b^2 + c^2)) (x^2 + y^2) +
b^2 ((b^2 + 2 (b^2 + c^2)) x^2 + (-3 b^2 +
2 (b^2 + c^2)) y^2)) x' y')/(x^2 + y^2)^(3/
2) + (a (c + Sqrt[b^2 + c^2]) y (c^4 (x^2 + y^2) +
3 c^3 Sqrt[b^2 + c^2] (x^2 + y^2) +
c Sqrt[b^2 + c^2] (7 b^2 + c^2) (x^2 + y^2) +
c^2 (2 b^2 + 3 (b^2 + c^2)) (x^2 + y^2) +
2 b^2 ((b^2 + 2 (b^2 + c^2)) x^2 - (b^2 -
2 (b^2 + c^2)) y^2)) (y')^2)/(x^2 + y^2)^(3/
2) + (8 a b^2 Sqrt[
b^2 + c^2] (c + Sqrt[b^2 + c^2])^2 y t' z')/
Sqrt[x^2 +
y^2] + (8 a b^2 (c + Sqrt[b^2 + c^2])^2 y^2 x' z')/(x^2 +
y^2)^(3/
2) - (8 a b^2 (c + Sqrt[b^2 + c^2])^2 x y y' z')/(x^2 +
y^2)^(3/2) - (4 a b^2 c (c + Sqrt[b^2 + c^2])^2 y (z')^2)/
Sqrt[x^2 + y^2]))};
To find the solution, I use the following procedure. First I take the values of the parameters to be:
b = 1
c = 20
a = 10^(-6)
Then I numerically evaluate them:
ngeodesic = N[geodesic];
And use the following procedure in order to obtain a solution for all the coordinates. First I calculate the first component of the four-velocity at parameter value equal to zero - this gives the remaining initial value for the NDSolve. After the NDSolve piece of code, I list the four-velocity vector's length to check if the system evolves the way it should.
minkowskimetric = {{-1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0,
0, 1}}
computeSoln[maxτi_, ivsi_, icsi_] :=
Block[{ivs, ics, i, χ, tmp, soln},
ics = icsi;
ivs = Join[{χ}, ivsi];
tmp = minkowskimetric;
tmp = ivs.(tmp.ivs);
χslv = Solve[tmp == uinvar, χ];
ivs[[1]] = Last[χ /. χslv];
deq = Table[
coords[[i]]''[τ] == (ngeodesic[[i]] /.
Join[Table[coords[[i]]' -> coords[[i]]'[τ], {i, 1, n}],
Table[coords[[i]] -> coords[[i]][τ], {i, 1, n}]]), {i, 1,
n}];
deq = Join[deq, Table[coords[[i]]'[0] == ivs[[i]], {i, 1, n}],
Table[coords[[i]][0] == ics[[i]], {i, 1, n}]];
soln = NDSolve[deq, coords, {τ, 0, maxτi}];
soln]
uinvar = -1;
udotu[solni_, τval_] := Block[{xα, uα},
xα =
Table[coords[[i]][τ] /. solni, {i, 1, n}] // Flatten;
uα = D[xα, τ];
xα = xα /. τ -> τval;
uα = uα /. τ -> τval;
uα.(minkowskimetric.uα)]
coordlist[τin_] :=
Table[ToString[coords[[i]]] <>
"=" <> {ToString[coords[[i]][τin] /. soln // First]}, {i, 1,
n}]
In what follows I use the initial values that actually yield some result. The behaviour looks like something obtained by my advisor, who used Maple. The ivs table is the initial velocities table - (x,y,z) components. The ics table is the table of initial coordinates (t,x,y,z).
maxτ = 40
ivs = {0.01, 0.01, 0.1};
ics = {0, 0.05, 0.05, 0.0};
soln = computeSoln[maxτ, ivs, ics];
ParametricPlot3D[
Evaluate[ {coords[[2]][τ], coords[[3]][τ],
coords[[4]][τ]}] /. soln, {τ, 0, maxτ}]
Plot[Evaluate[Table[coords[[i]][τ] /. soln, {i, 1, n}]], {τ,
0, maxτ}, PlotStyle -> Table[Hue[(i - 1)/n], {i, 1, n}],
PlotLegends -> coords, AxesLabel -> {"τ", "Coordinate"}]
ParametricPlot[
Evaluate[ {coords[[2]][τ], coords[[3]][τ]}] /.
soln, {τ, 0, maxτ}]
Join[{{"", "", "", "u*u values"}},
Table[{"τ=", ToString[i], "->", udotu[soln, i]}, {i, 0,
maxτ, maxτ/5}]] // TableForm
The problem begins when I try to change the initial conditions to, let's say:
ics = {0, 0.5, 0.05, 0.0};
In that case (and literally any other) I get the error:
NDSolve::ndsz: At τ == 0.34526573966566704`, step size is effectively zero; singularity or stiff system suspected.
My advisor (who does not use Mathematica) gets very nice results for almost any imaginable parameters - since the physical problems only occur on the z coordinate axis. I am at loss at the moment - I suspect that the very small values (e.x. the amplitude "a") are the cause for the problem. However, I have no idea if I could even rescale all the quantities used here - the problem is, after all, highly nonlinear, even in the parameters a,b,c.
EDIT:
I put the whole code here: besselgravcart
NDSolvebecomes numerically unstable forics[[2]]larger than about0.1189.Method->"StiffnessSwitching"does not help. Neither does increasingWorkingPrecision(with all inputs rationalized). Perhaps, better choices of other options would. By the way, comments by other experienced users suggest that Maple is better than Mathematica at solving some ODEs. $\endgroup${coords[[2]][τ], coords[[3]][τ]}space? Alternatively, does the solution satisfy any conservation laws? $\endgroup$x - yspace come in from larger, skim a circle of constantrfor some distance and then go out to largeragain. So, I speculate that this circle of constantris a separatrix. If so, then infinitesimal errors in the computation near the separatrix would cause the computation to fail catastrophically. First integrals of the motion might help to alleviate this problem. Can the equations be rewritten conveniently in terms orrandthetainstead ofxandy(i.e., in cylindrical coordinates)? In other words, is the metric axisymmetric? $\endgroup$tis a linear function ofz, so one of the four ODEs is not needed. Correct? Is the metric also axisymmetric, which could reduce another equation to first order or eliminate it entirely? By the way, how many cycles inx - ydo you need to compute to obtain a useful answer? $\endgroup$