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

Question

我读了费曼物理讲义第9章并尝试了自己的模拟。我用黎曼积分来计算速度和位置。虽然所有的起点都是一样的，但我的轨道看起来像一条双曲线。这是讲义：https://www.feynmanlectures.caltech.edu/I_09.html (Table 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 越大，"local linearization error" 就越大（可能有一个更好的数学术语...）。

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说明你在做什么：

A 'brute force' kind of numerical integration can be done, if the integrand is reasonably well-behaved (i.e. piecewise continuous and of bounded variation), by evaluating the integrand with very small increments.

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

An important part of the analysis of any numerical integration method is to study the behavior of the approximation error as a function of the number of integrand evaluations.

有many smart approaches to make better guesses and use larger time steps. You might have heared of the Runge–Kutta method to solve differential equations. It seems to become Simpson's rule mentioned in the link above for non-differential equations, see here.

Answer 2

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

这个问题可以有两种解法：

将时间间隔减小到更小的值，这样即使在整个过程中错误放大后也不会与实际解决方案有很大偏差。现在要注意的是，如果您减少时间间隔 (dt) 而不是增加步数，那么您可以在更短的时间内看到系统的演变。因此，随着时间间隔 (dt) 的减少，您还必须增加步数。 我检查了你的代码，似乎如果你将 dt = 0.0001 并将步数设为 100000 而不是 1000，你就会得到你正在寻找的美丽椭圆轨道。 另外，删除或注释掉 plt.xlim 和 plt.ylim 行以清楚地查看您的情节。
实现微分方程数值解的龙格-库塔方法。这种方法每次迭代都可以更好地收敛到实际解决方案。因为，这将花费更多时间并更改您的代码，这就是为什么我建议将其作为第二个选项，否则此方法比 Euler 的方法更高级和更通用。

注意： 即使没有任何更改，解也不会表现得像双曲线。对于您为系统提供的初始值，解决方案给出了一条有界曲线，但仅仅因为误差放大，它正在螺旋上升到一个点。如果您将步数增加到 10000 并设置 dt = 0.01.

，您会注意到这种螺旋式上升

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

Numerical integration: Why does my orbit simulation yield the wrong result?

python

physics

numerical-methods

astronomy

physics-engine