lagrange-poisson 0.2.0

拉格朗日-泊松问题数值求解器 / Lagrange-Poisson Problem Numerical Solver

拉格朗日-泊松问题求解器

一个用于求解刚体绕定点在重力场中运动的Python数值计算工具包，提供图形界面和多种数值方法。

功能特点

• 实现三种经典数值方法：

  • 显式欧拉法 (一阶精度)
  • 改进欧拉法/Heun法 (二阶精度)
  • 四阶Runge-Kutta法 (四阶精度)

• 图形用户界面 (GUI)

  • 基于 Tkinter 的直观界面
  • 实时参数调整
  • 结果可视化展示
  • 数据导出功能

• 科学可视化

  • 基于 Matplotlib 的高质量绘图
  • 角速度和方向余弦的时间演化曲线

• 完整的物理模型

  • 欧拉动力学方程
  • 泊松运动学方程
  • 支持任意主转动惯量配置
  • 支持任意质心位置

• 欧拉角输出

  • 泊松变量到欧拉角的自动转换
  • 支持进动角φ、章动角θ、自转角ψ的输出
  • GUI中可选显示模式（泊松变量/欧拉角）

安装

pip install lagrange-poison

或从源码安装：

git clone https://github.com/yourusername/lagrange-poison.git
cd lagrange-poison
pip install -e .

快速开始

使用命令行启动GUI

lagrange-poison

在Python代码中使用

import numpy as np
from lagrange_poison import LagrangePoissonSolver, RK4Method

# 创建求解器
solver = LagrangePoissonSolver(
    A=1.0,   # 主转动惯量
    B=2.0,
    C=3.0,
    x0=0.1,  # 质心坐标
    y0=0.0,
    z0=0.0,
    m=1.0,   # 质量
    g=9.81   # 重力加速度
)

# 设置初始条件
initial_conditions = np.array([0.5, 0.1, 0.05, 0.0, 0.0, 1.0])
# [p0, q0, r0, gamma1_0, gamma2_0, gamma3_0]

# 求解
method = RK4Method()
t, y = solver.solve(
    initial_conditions=initial_conditions,
    t_span=(0, 20),
    dt=0.01,
    method=method
)

# 绘制结果
import matplotlib.pyplot as plt

fig, axes = plt.subplots(2, 3, figsize=(12, 8))
labels = ['p', 'q', 'r', 'γ₁', 'γ₂', 'γ₃']
titles = ['角速度 p', '角速度 q', '角速度 r',
          '方向余弦 γ₁', '方向余弦 γ₂', '方向余弦 γ₃']

for i, ax in enumerate(axes.flat):
    ax.plot(t, y[:, i])
    ax.set_title(titles[i])
    ax.set_xlabel('时间 (s)')
    ax.set_ylabel(labels[i])
    ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

泊松变量与欧拉角转换

本软件支持将泊松变量 $(\gamma_1, \gamma_2, \gamma_3)$ 转换为欧拉角 $(\phi, \theta, \psi)$：

# 转换为欧拉角
euler = solver.poisson_to_euler_angles(t, y, phi0=0.0)

phi = euler[:, 0]   # 进动角
theta = euler[:, 1] # 章动角
psi = euler[:, 2]   # 自转角

# 绘制欧拉角
fig, axes = plt.subplots(2, 3, figsize=(12, 8))
labels = ['φ (rad)', 'θ (rad)', 'ψ (rad)', 'φ (°)', 'θ (°)', 'ψ (°)']
titles = ['进动角', '章动角', '自转角'] * 2

for i in range(3):
    # 弧度
    axes[0, i].plot(t, euler[:, i])
    axes[0, i].set_title(titles[i])
    axes[0, i].set_ylabel(labels[i])

    # 角度
    axes[1, i].plot(t, np.degrees(euler[:, i]))
    axes[1, i].set_title(titles[i])
    axes[1, i].set_ylabel(labels[i + 3])
    axes[1, i].set_xlabel('时间 (s)')

plt.tight_layout()
plt.show()

转换公式：

• 章动角: $\theta = \arccos(\gamma_3)$
• 自转角: $\psi = \operatorname{atan2}(\gamma_1, \gamma_2)$
• 进动角: $\phi_{n+1} = \phi_n + h\frac{p_n\gamma_{1,n} + q_n\gamma_{2,n}}{\gamma_{1,n}^2 + \gamma_{2,n}^2}$

注意： 欧拉角表示在 $\sin\theta \approx 0$ 时存在奇异点，这是欧拉角表示法的固有特性。

物理背景

本软件求解的是拉格朗日-泊松问题，即刚体绕定点在重力场中的运动。

数学方程

欧拉动力学方程：

Aṗ + (C-B)qr = M₁
Bq̇ + (A-C)rp = M₂
Cṙ + (B-A)pq = M₃

泊松运动学方程：

γ̇₁ = rγ₂ - qγ₃
γ̇₂ = pγ₃ - rγ₁
γ̇₃ = qγ₁ - pγ₂

其中：

• p, q, r 是角速度在体轴系中的分量
• γ₁, γ₂, γ₃ 是重力方向在体轴系中的方向余弦
• A, B, C 是主转动惯量
• M₁, M₂, M₃ 是重力矩分量

参数说明

参数	说明	单位
A, B, C	主转动惯量	kg·m²
x₀, y₀, z₀	质心在体坐标系中的位置	m
m	刚体质量	kg
g	重力加速度	m/s²
p₀, q₀, r₀	初始角速度	rad/s
γ₁₀, γ₂₀, γ₃₀	初始方向余弦	无量纲

数值方法

1. 显式欧拉法

最简单的数值方法，一阶精度：

y_{n+1} = y_n + h·f(y_n)

2. 改进欧拉法 (Heun方法)

二阶精度的预测-校正方法：

k₁ = f(y_n)
k₂ = f(y_n + h·k₁)
y_{n+1} = y_n + (h/2)(k₁ + k₂)

3. 四阶Runge-Kutta法 (RK4)

最常用的四阶精度方法：

k₁ = f(y_n)
k₂ = f(y_n + h·k₁/2)
k₃ = f(y_n + h·k₂/2)
k₄ = f(y_n + h·k₃)
y_{n+1} = y_n + (h/6)(k₁ + 2k₂ + 2k₃ + k₄)

依赖项

• Python >= 3.8
• NumPy >= 1.20.0
• Matplotlib >= 3.3.0

开发

运行测试

pip install -e ".[dev]"
pytest

构建和发布

# 构建
python -m build

# 发布到PyPI (需要配置 API token)
python -m twine upload dist/*

许可证

MIT License

作者

Your Name [e-mail]

参考文献

1. Arnold, V. I. (1989). Mathematical Methods of Classical Mechanics
2. Landau, L. D., & Lifshitz, E. M. (1976). Mechanics
3. Goldstein, H. (1980). Classical Mechanics

贡献

欢迎提交 Issue 和 Pull Request！