数值积分:为什么我的轨道模拟会产生错误的结果? [英] Numerical integration: Why does my orbit simulation yield the wrong result?

问题描述

我阅读了费曼物理讲座第 9 章，并尝试了我自己的模拟.我使用黎曼积分来计算速度和位置.虽然所有的起点都一样，但我的轨道看起来像双曲线.这是讲义:https://www.feynmanlectures.caltech.edu/I_09.html(表 9.2)

I read Feynman's Lecture on Physics Chapter 9 and tried to my own simulation. I used Riemann integrals to calculate velocity and position. Although all start-entry is same, my orbit look's like a hyperbola. Here is lecture note: 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()

推荐答案

dt 应该无限小，以便集成工作.

dt is supposed to be infinitly small for the integration to work.

dt 越大，局部线性化误差"就越大(可能有一个更好的数学术语......).

The bigger dt the bigger the "local linearization error" (there is probably a nicer mathematical term for that...).

0.1 for dt 可能看起来足够小，但是为了让您的结果正确收敛(朝向现实)，您还需要检查更小的时间步长.如果较小的时间步长收敛于相同的解，您的方程足够线性以使用更大的时间步长(并节省计算时间.)

0.1 for dt may look small enough, but for your result to converge correctly (towards reality) you need to check smaller time steps, too. If smaller time steps converge towards the same solution your equation is linar enough to use a bigger time step (and save comptation time.)

用

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

和

for i in range(0, 100000):
    earth.calculateOrbit(0.001, planets)

在两次计算中，自开始以来经过的总时间与您的原始值相同.但结果却大不相同.因此，您可能必须使用更小的 dt.

In both calculations the overall time that has passed since the beginning is the same as with your original values. But the result is quite different. So you might have to use an even smaller dt.

更多信息:

https://en.wikipedia.org/wiki/Trapezoidal_rule

并且此页面说明您在做什么:

可以进行蛮力"数值积分，如果被积函数的行为相当好(即分段连续和有界变化)，通过评估被积函数非常小增量.

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.

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

And what I tried to explain above:

任何数值积分方法分析的重要组成部分是研究作为函数的近似误差的行为被积函数评估的数量.

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.

有许多智能方法来做出更好的猜测并使用更大的时间步长.您可能听说过用于求解微分方程的 Runge–Kutta 方法.对于非微分方程，它似乎成为上面链接中提到的 Simpson's rule，参见 此处.

There are 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.

这篇关于数值积分:为什么我的轨道模拟会产生错误的结果?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！