我想数值求解维基百科在 Matlab 中给出的克尔测地线的轨迹方程。

轨迹如下所示：

我实现了方程并用标准的龙格-库塔方法求解。这是我的代码：

% Initial parameter (r, theta, phi, t)
r       = x(1);
theta   = x(2);
phi     = x(3);
t       = x(4);

% Constants of motion
mu = const(1);
E = const(2);
L = const(3);
Q = const(4);

% Abbreviations
Sigma = r^2 + a^2 * cos(theta)^2;
Delta = r^2 - 2 * M * r + a^2;
K = Q + (L - a*E)^2;
X = a^2*(E^2 + mu) - L^2 - Q;
R = (E^2 + mu)*r^4 - 2*M*mu*r^3 + X*r^2 + 2*M*K*r -a^2*Q;
O = Q + cos(theta)^2*(a^2*(E^2 + mu) - L^2/sin(theta)^2);

% ODE from O Neill Chapter 4.2.2 Geodesics
dx(1) = sqrt(R / Sigma^2);
dx(2) = sqrt(O / Sigma^2);
dx(3) = ((L - a*sin(theta)^2*E) / sin(theta)^2 + a*((r^2 + a^2)*E - L*a)/(Delta) ) / Sigma^2;
dx(4) = (a*(L - a*E*sin(theta)^2) + (r^2 + a^2)*((r^2 + a^2)*E - L*a)/(Delta)) / Sigma^2;

我对 R 项和 Q 项有点困惑，因为有一个平方根。我如何在计算中处理这些符号？

我知道存在具有动量分量的方程的另一种形式，但我不想使用这种方法。