计算科学 - 为初值 ODE 实现欧拉方法 - 吾爱随笔录

为初值 ODE 实现欧拉方法

计算科学软件颂

2021-12-13 01:51:24

在我的物理课上，我必须计算以，质量）离开表面的角度（半径R ，质量 M ）发射（非常快）的射弹的轨迹。弹丸将达到最大距离并开始围绕地球旋转。问题是它是否会撞击地球。我们通过分析解决了这个问题，但我想知道如何以数字方式模拟这样的问题。 $v_0$ $R$ $M$

所以基本上我从

{\vec{r}}_{0} = (\begin{matrix} 0 \\ 0 \\ R \end{matrix})

$\vec{r}_0 = \begin{pmatrix} 0 \\ 0 \\ R \end{pmatrix}$

和

{\vec{v}}_{0} = v_{0} (\begin{matrix} \cos (30 °) \\ 0 \\ \sin (30 °) \end{matrix})

$\vec{v}_0 = v_0 \begin{pmatrix} \cos(30°) \\ 0 \\ \sin(30°) \end{pmatrix}$

弹丸的加速度只是

\vec{a} (\vec{r}) = \frac{\vec{F} (\vec{r})}{m} = G \frac{M}{r^{3}} \vec{r}

$\vec{a}(\vec{r}) = \frac{\vec{F}(\vec{r})}{m} = G \frac{M}{r^3} \vec{r}$

我可以实现一个简单的 C++/Python 程序，它基本上可以：

{\vec{r}}_{t + 1} = {\vec{r}}_{t} + {\vec{v}}_{t} \cdot d t

$\vec{r}_{t+1} = \vec{r}_{t} + \vec{v}_{t} \cdot \mathrm dt$

{\vec{v}}_{t + 1} = {\vec{v}}_{t} + {\vec{a}}_{t} \cdot d t

$\vec{v}_{t+1} = \vec{v}_{t} + \vec{a}_{t} \cdot \mathrm dt$

但这似乎有点不适合，因为在这两种语言中我都没有“自然”向量类。

Octave 或 Mathematica 会更适合这项任务吗？为了在 gnuplot 中绘制轨迹，我将如何实现这个问题以获得 x、y、z 坐标列表？

3个回答

只需使用 Python + NumPy。这是代码

from math import pi, sin, cos
from numpy import array, vstack
from numpy.linalg import norm

from pylab import plot, savefig

R = 10.
M = 1.
G = 1.
phi = 30 * 180. / pi
v0_magnitude = 0.1
r = array([0, 0, R])
v = v0_magnitude * array([cos(pi), 0, sin(pi)])
trajectory = []

dt = 0.1
t = 0
while t < 100:
    trajectory.append(r.copy())
    a = -G * M * r / norm(r)**3
    v += a * dt
    r += v * dt
    t += dt

trajectory = vstack(trajectory)
plot(trajectory[:, 0], trajectory[:, 2])
savefig("trajectory.png")

这是图表：

轨迹.png

我还在这个gist上传了上面的代码和图表。

对于涉及随时间演变的系统，一种选择是考虑使用 MATLAB 或 Modelica。MATLAB 当然是专有的，但 Modelica 存在于许多实现中，其中一些，例如OpenModelica，绝对是开源的。

作为一个对 Modelica 比较陌生的人（我大约一年前才知道它），我发现它非常擅长实现基于动力系统的问题：常微分方程组或微分代数方程组。特别是，问题的设置是高度模块化的，您可以以“自然”形式编写方程，而不必以任意形式重写它们以将它们按摩成语言或实现所需的任何格式，例如 MATLAB 要求的。OpenModelica 还提供了一些集成的绘图功能。

Maple 和 Mathematica 都有简单的接口来数值求解微分方程并绘制结果。这里 http://www.math.tamu.edu/~bangerth/teaching/2010-fall-442/project-tools.pdf http://www.math.tamu.edu/~bangerth/teaching/2010-fall- 442/project-tools.mw 是 Maple 的一个示例，带有很酷的图表。

其它你可能感兴趣的问题

上一篇当候选点在凸包上时，是否有一种特殊的算法来计算凸包排序？下一篇并行求解 PDE