这段代码
k4r = f1 (t(我)+ h r (i) + k3r * h, rdot (i) + k3rdot * h,θ(i) + k3theta * h, thetad(我)+ k3r * h);
k4rdot = f2 (t(我)+ h r (i) + k3r * h, rdot (i) + k3rdot * h,θ(i) + k3theta * h, thetad(我)+ k3r * h);
k4theta = f3 (t(我)+ h r (i) + k3r * h, rdot (i) + k3rdot * h,θ(i) + k3theta * h, thetad(我)+ k3r * h);
k4thetad = f4 (t(我)+ h r (i) + k3r * h, rdot (i) + k3rdot * h,θ(i) + k3theta * h, thetad(我)+ k3r * h);
应该是这个
k4r = f1 (t(我)+ h r (i) + k3r * h, rdot (i) + k3rdot * h,θ(i) + k3theta * h, thetad(我)+ k3thetad * h);
K4rdot = f2(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3thetad*h);
k4 = f3(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3thetad*h);
K4thetad = f4(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3thetad*h);
话虽如此，你的代码是
道路
它太复杂了。你基本上需要为四个标量方程写四组RK4代码。它会是
多
最好是定义一个4元素的状态向量，然后写
一个
RK4方程的集合来处理传播。这意味着你目前所写的代码的1/4，更少的机会出现像你所做的复制粘贴错误。
还有这些初始条件
R (1) = 1000;卫星的初始位置
Rdot (1) = 50;卫星的初始转译率
(1) = 0;卫星的初始角位置
Thetad (1) = 0;卫星初始角速率
首先把卫星放在地球表面内，结果是直线运动。我觉得这不是你想要的。
此外,
Q1:为什么f2 - g *Me/r^2 - r*thetad^2不是
编辑:
错误符号修正建议如上。我应该建议的是保持f2不变，并进行更改
为什么不是f4 -2*rdot*thetad/r
另外，你可以看到这个相关的线程:
