使用scipy odeint求解二体问题时遇缺失位置参数错误求助

问题描述 #

尝试用scipy的odeint求解等质量引力二体问题，编写的代码如下：

#N Body test case

#%%

#import modules

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

#%%

#define constants

G=6.67e-11
AU=1.496e11
m1= m2=1.989e30

def twobody(x1, y1, vx1, vy1, x2, y2, vx2, vy2):
    
    rx=x1-x2
    ry=y1-y2
    r=np.sqrt(rx**2+ry**2)
    
    f=[vx1, vy1, vx2, vy2, G*m2*rx/r**3, G*m2*ry/r**3, G*m1*rx/r**3, G*m2*ry/r**3]
    
    return f

t=np.linspace(0, 1.577e8, 2000)

initial=[-0.5*AU, -0.5*AU, 0,0,0,0,-15000, 15000]

twobodysol=odeint(twobody, initial, t, args=(m1, m2))

运行后收到错误：

twobody() missing 4 required positional arguments: 'x2', 'y2', 'vx2', and 'vy2'

错误原因分析 #

1. odeint函数参数规范不匹配：odeint要求被积分的函数必须遵循func(y, t, *args)的签名，其中y是包含所有状态变量的一维数组，t是当前时间点，*args是额外常量参数。你的twobody函数把状态变量拆分成多个独立参数，odeint无法正确传递参数，导致报错。
2. 未接收额外参数：你在odeint的args中传入了m1,m2，但twobody函数的参数列表里没有这两个变量，即使参数规范正确也会报错。
3. 加速度方向错误：引力是吸引力，当前代码中计算的加速度方向不符合物理规律，会导致后续轨道计算错误。

修正后的代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 定义常量
G = 6.67e-11
AU = 1.496e11
m1 = m2 = 1.989e30

def twobody(y, t, G, m1, m2):
    # 从状态向量中拆分所有变量
    x1, y1, vx1, vy1, x2, y2, vx2, vy2 = y
    
    # 计算相对位置与距离
    rx = x1 - x2
    ry = y1 - y2
    r = np.sqrt(rx**2 + ry**2)
    
    # 计算加速度：引力为吸引力，方向指向对方
    ax1 = -G * m2 * rx / r**3  # m1的x方向加速度，指向m2
    ay1 = -G * m2 * ry / r**3  # m1的y方向加速度
    ax2 = G * m1 * rx / r**3   # m2的x方向加速度，指向m1
    ay2 = G * m1 * ry / r**3   # m2的y方向加速度
    
    # 返回状态向量的导数：位置的导数是速度，速度的导数是加速度
    return [vx1, vy1, ax1, ay1, vx2, vy2, ax2, ay2]

# 时间数组：覆盖约5年（1.577e8秒），取2000个采样点
t = np.linspace(0, 1.577e8, 2000)

# 初始状态：[x1, y1, vx1, vy1, x2, y2, vx2, vy2]，调整为对称形式便于观察轨道
initial = [-0.5*AU, 0, 0, 15000, 0.5*AU, 0, 0, -15000]

# 调用odeint求解：args参数顺序需与twobody的额外参数一致
twobodysol = odeint(twobody, initial, t, args=(G, m1, m2))

# 绘制轨道
plt.figure(figsize=(8,8))
plt.plot(twobodysol[:,0]/AU, twobodysol[:,1]/AU, label='m1')
plt.plot(twobodysol[:,4]/AU, twobodysol[:,5]/AU, label='m2')
plt.scatter(0,0, color='yellow', label='质心')  # 等质量二体质心在原点
plt.xlabel('X (AU)')
plt.ylabel('Y (AU)')
plt.legend()
plt.axis('equal')
plt.title('等质量二体问题轨道')
plt.show()

修正说明 #

• 调整twobody函数的参数结构，符合odeint的调用规范，从状态向量中拆分出各个变量。
• 在函数参数中添加G,m1,m2，并在odeint的args中按顺序传入这些常量。
• 修正加速度的符号，确保引力方向正确（两个物体相互吸引）。
• 调整初始状态为对称形式，更直观地观察等质量二体的绕质心运动。

内容的提问来源于stack exchange，提问作者XBlake97