数字积分：为什么我的轨道模拟产生错误的结果？

Question

我阅读了费曼的物理讲座第9章，并尝试进行自己的模拟。我使用黎曼积分来计算速度和位置。尽管所有开始输入都是相同的，但我的轨道看上去像个双曲线。这是讲义：https://www.feynmanlectures.caltech.edu/I_09.html（表9.2）

import time
import matplotlib.pyplot as plt
x=list()
y=list()
x_in=0.5
y_in=0.0
x.append(x_in)
y.append(y_in)

class Planet:
    def __init__(self,m,rx,ry,vx,vy,G=1):
        self.m=m
        self.rx=rx
        self.ry=ry
        self.a=0
        self.ax=0
        self.ay=0
        self.vx=vx
        self.vy=vy
        self.r=(rx**2+ry**2)**(1/2)
        self.f=0
        self.G=1
        print(self.r)
    def calculateOrbit(self,dt,planet):
        self.rx=self.rx+self.vx*dt
        self.vx=self.vx+self.ax*dt
        self.ax=0
        for i in planet:
            r=((((self.rx-i.rx)**2)+((self.ry-i.ry)**2))**(1/2))
            self.ax+=-self.G*i.m*self.rx/(r**3)

        self.ry=self.ry+self.vy*dt
        self.vy=self.vy+self.ay*dt
        self.ay=0
        for i in planet:
                self.ay+=-self.G*i.m*self.ry/(r**3)
        global x,y
        x.append(self.rx)
        y.append(self.ry)

        #self.showOrbit()
    def showOrbit(self):
        print(""" X: {} Y: {} Ax: {} Ay: {}, Vx: {}, Vy: {}""".format(self.rx,self.ry,self.ax,self.ay,self.vx,self.vy))

planets=list();
earth = Planet(1,x_in,y_in,0,1.630)
sun =   Planet(1,0,0,0,0)
planets.append(sun)
for i in range(0,1000):
    earth.calculateOrbit(0.1,planets)

plt.plot(x,y)
plt.grid()
plt.xlim(-20.0,20.0)
plt.ylim(-20.0,20.0)
plt.show()

Answer 1

dt应该很小，才能使集成正常工作。

dt越大，“局部线性化误差”越大（可能存在一个更好的数学术语……）。

0.1等于dt可能看起来很小，但是为了使结果正确收敛（朝着现实），您还需要检查较小的时间步长。如果较小的时间步长收敛于相同的解，则您的方程式为linar enough，以使用较大的时间步长（并节省计算时间。）

尝试使用]进行代码>

for i in range(0, 10000):
    earth.calculateOrbit(0.01, planets)
和

for i in range(0, 100000):
    earth.calculateOrbit(0.001, planets)
在两次计算中，从开始算起经过的总时间与原始值相同。但是结果却大不相同。因此，您可能必须使用更小的dt。

更多信息：

https://en.wikipedia.org/wiki/Trapezoidal_rule
[this page说明您在做什么：

如果 被积数合理地表现良好（即分段连续且具有 有界变异），通过对被积物进行非常小的评估 增量。

以及我上面试图解释的内容：

任何数值积分方法分析的重要部分 是研究近似误差的行为与 被评估数。

many smart approaches可以做出更好的猜测并使用较大的时间步长。您可能听说过Runge–Kutta method用于求解微分方程。对于非微分方程，它似乎已成为上面链接中提到的Simpson's rule，请参阅here。

Answer 2

问题是微分方程的数值积分方法或数值解。您正在使用的方法（用于微分方程数值解的欧拉方法），尽管它与实际值非常接近，但误差仍然很小。如果在多次迭代中都使用了这个略有错误的值（例如您已完成1000个步骤），则此错误在每一步都会变得越来越大，从而给您错误的结果。

可以有两个解决方案两个这个问题：