在Mathematica中模拟星体运动

2020-05-02

引言

我们想看一看,星体互相之间受到吸引,那么这些星体的运动轨迹是怎样.我们知道,牛顿第二定律指出:物体的加速度与其受到的合力同方向,且加速度的大小与合力的大小成正比,与自身的质量成反比,而万有引力定律指出:两个物体总是互相吸引,或者说两个物体之间存在某种引力,该引力的方向平行于两个物体的连线,引力的大小正比于两物体的质量的乘积,反比于两物体的质心距离的平.模拟星体的运动,我们主要就是根据这两条物理定律.

要用到的基本定理的形式化

我们打算在计算机中模拟星体的运动,实际上,就是模拟星体随着时间变化不断地受力-加速-移动的这个过程,其中用到的牛顿第二定律,我们记加速度为$\overrightarrow{a}$,记质量为$m$,记物体所受到的力的和(合力)为$\overrightarrow{F}$,写起来就是:

$$ \overrightarrow{a} = \frac{1}{m} \overrightarrow{F} $$

而另外一条同样重要的要用到的物理定律是万有引力定律,树上的苹果掉到地上受到的是万有引力(更常见的叫做重力),地球围绕着太阳转,也是因为太阳对地球的万有引力在一直「拉着」地球,使得地球不至于因为自身惯性飞太阳系外边去,那么,对于万有引力,我们记第一个物体的质量为$m_1$,第二个物体的质量为$m_2$,万有引力常量为$G$,两物体的质心的距离是$r$,那么,两物体之间的万有引力的大小为:

$$ F = \frac{G m_1 m_2}{r^2} $$

至于方向呢,可以看做是物体1拉着物体2,这时引力的方向是从物体2指向物体1的,作用在物体2身上,也可以看做是物体2拉着物体1,这时引力的方向是从物体1指向物体2的,作用在物体1身.实际计算中,我们选择其中一个物体作为「中心」,另外一个物体作为「边沿」物体,中心的地位类似于荷包蛋中间的蛋黄,我们只需把万有引力看做是中心物体拉着边沿物体的力即可,并不需要先后计算物体1对物体2的作用,再计算物体2对物体1的作用,因为,物体1拉着物体2移动了一段距离,也相当于物体2拉着物体1移动了一段距离.

基本的变速运动

在开始行星运动轨迹的计算之前,我们来看一个更简单的同时也是更基本的问题:一个物体在$t_0$时刻的位置是在$\overrightarrow{x_0}$,同时,该物体在$t_0$时刻的速度是$\overrightarrow{v_0}$,假设经过了非常短的时间$\Delta t$,使得在$t_1 = t_0 + \Delta t$时刻,物体的速度变成了$\overrightarrow{v_1}$,试估计该物体在$t_0$到$t_1$期间行进的轨迹.

估算物体的运动轨迹

估算物体的运动轨迹

一种朴实的解决该问题的思路会考虑使用微元法,其实也就是,把物体在$t_0$到$t_1$期间的运动,看成是非常非常多个匀速直线运动的组合,更加具体地,我们会这样来考虑,我们把$t_0$到$t_1$的时间间隔$\Delta t$分成$n$份,那么,在$t_0$到$t_0 + \Delta t/n$,物体的运动可近似看做是匀速直线运动,而在$t_0 + \Delta t/n$到$t_0 + 2 \Delta t/n$,物体的运动也同样可看做是匀速直线运动,就这样,一直到$t_0 + (n-1) \Delta t/n$和$t_0 + n \Delta t/n$也就是$t_1$这个区间,物体的在这里列出的每一个区间都是匀速直线运.同时,我们还把从时刻$t_0$到时刻$t_1$速度的变化量$\Delta \overrightarrow{v} = \overrightarrow{v_1} - \overrightarrow{v_0}$也分成$n$份,那么,在$t_0$到$t_0 + \Delta t/n$期间,物体的运动速度近似等于$\overrightarrow{v_0}$,在$t_0 + \Delta t /n$到$t_0 + 2 \Delta t /n$期间,物体的运动速度近似等于$\overrightarrow{v_0} + \Delta \overrightarrow{v} / n$,在$t_0 + 2 \Delta t / n$期间,物体的运动速度近似等于$\overrightarrow{v_0} + 2 \Delta \overrightarrow{v} / n$,就这样,由于总的时间间隔$\Delta t = t_1 - t_0$是非常小的,所以,可以用这种等差数列来估计每一个长度为$\Delta t / n$的子区间的速.其实也就是,对于这个$t_0$到$t_1$的运动轨迹问题,我们把$\Delta t$分成$n$份,然后一段一段地求出这个轨迹.

更加简洁地,我们可以把该物体从$t_0$时刻到$t_1$时刻的运动,看成是$n$次物理状态的迭代,把该物体在$t_0$时刻的物理状态记做$(\overrightarrow{x^{(0)}}, \overrightarrow{v^{(0)}})^{\mathsf{T}}$,其中的$\overrightarrow{x}$表示位置,$\overrightarrow{v}$表示速度,而上标$(0)$则表示是第几次迭代,$\mathsf{T}$代表矩阵的转置,那么,物体的运动,实际上可以看做是物理状态的不断迭代:

$$ \begin{bmatrix} \overrightarrow{x^{(0)}} \\
\overrightarrow{v^{(0)}} \end{bmatrix} \rightarrow \begin{bmatrix} \overrightarrow{x^{(1)}} \\
\overrightarrow{v^{(1)}} \end{bmatrix} \rightarrow \begin{bmatrix} \overrightarrow{x^{(2)}} \\
\overrightarrow{v^{(2)}} \end{bmatrix} \rightarrow \cdots \rightarrow \begin{bmatrix} \overrightarrow{x^{(n)}} \\
\overrightarrow{v^{(n)}} \end{bmatrix} $$

具体地,这个迭代的过程是这样的:

$$ \begin{cases} \overrightarrow{x^{(n)}} = \overrightarrow{x^{(n-1)}} + \Delta t_{\epsilon} \overrightarrow{v^{(n-1)}} \\
\overrightarrow{v^{(n)}} = \overrightarrow{v^{(n-1)}} + \Delta \overrightarrow{v} / n \end{cases} $$

式中,$\Delta \overrightarrow{v} = \overrightarrow{v_1} - \overrightarrow{v_0}$,而$\Delta t_{\epsilon} = \Delta t / n = (t_1 - t_0) / n$,并且$\overrightarrow{x^{(0)}} = \overrightarrow{x_0}$,$\overrightarrow{v^{(0)}} = \overrightarrow{v_0}$,那么,显然,物体在$t_0$时刻到$t_1$时刻期间的运动轨迹,其实就可以形式化地表示为

$$ \overrightarrow{x^{(0)}} \rightarrow \overrightarrow{x^{(1)}} \rightarrow \overrightarrow{x^{(2)}} \rightarrow \cdots \rightarrow \overrightarrow{x^{(n)}} $$

上列式子表示,物体先在$\overrightarrow{x^{(0)}}$的位置,然后在$\overrightarrow{x^{(1)}}$的位置,然后在$\overrightarrow{x^{(2)}}$的位置,……,最后在$\overrightarrow{x^{(n)}}$的位置,而迭代方程,则告诉了我们这些个$\overrightarrow{x^{(i)}}, i = 0,1,2,\cdots,n$具体是怎么来的.

符号化的轨迹表示

符号化的轨迹表示

最终,我们只需按照坐标把$\overrightarrow{x^{(i)}},i=1,\cdots,n$一个个的描出来,就得到物体的运动轨迹啦.

接下来,我们尝试用Mathematica语言在Mathematica软件中实现和验证我们的推理:

move[position0_, velocity0_, velocity1_, deltaT_] := (
    noCuts = 10;
    deltaV = velocity1 - velocity0;
    smallMove[{oldPosition_, oldVelocity_}] := (
        newPosition = oldPosition + oldVelocity*((deltaT)/(noCuts));
        newVelocity = oldVelocity + (deltaV)/(noCuts);
        {newPosition, newVelocity}
    );
    #[[1]] & /@ NestList[smallMove, {position0, velocity0}, noCuts]
);

实现了这样的一个函数,它接受一个初始位置,初速度和末速度,以及运动时间,经过计算之后输出物体的运动轨.先来看一下,假如说一个物体初始时刻的坐标是$(1,2)$,初速度和末速度都是$(2,1)$,经过了长度为$4$的时间,求物体的运动轨.也就是:

画出运动轨迹

画出运动轨迹

这其实是一个匀速直线运动,按照我们学过的方法,物体在最终时刻的位置等于

$$ (1, 2) + 4 \times (2, 1) = (9, 6) $$

可以看到程序的输出和我们的计算是符合的,并且,我们还用作图函数画出了物体的运动轨迹.

下面,我们来看另外一个例子:物体的初始位置仍然是在$(1,2)$,在长度为$0.5$的时间内,物体的速度由$(3,1)$变为$(2,3)$,求物体的运动轨.用我们实现了的函数来求解就是:

轨迹作图

轨迹作图

我们写的程序给出的结果是正确的吗?为了验证,我们可以分别验证点在$x$轴的投影的运动轨迹和点在$y$轴的投影的运动轨迹,我们学过,假设一个物体的速度从$v_0$变为$v_1$,如果这个加速过程的加速度大小恒定在$a$,该物体在这段加速过程中行进的距离记做$x$,那么有

$$ v_1^2 - v_0^2 = 2ax $$

或者,当该式子用来比较理论和实验观察值时

$$ |v_1^2 - v_0^2 - 2ax| \leq \epsilon $$

其中的$\epsilon$表示误差,一般是一个非常小的正.那么对于物体在$x$轴的投影的运动,有

$$ |2^2 - 3^2 - 2 \times ((2-3)/0.5) \times (2.2625-1)| = 0.05 $$

而对于物体在$y$轴的投影的运动,有

$$ |3^2 - 1^2 - 2 \times ((3-1)/0.5) \times (2.975-2)| = 0.2 $$

观察到有一定的误差,接下来我们试着增大区间的个数,然后再重新算一遍:

move[position0_, velocity0_, velocity1_, deltaT_] := (
    noCuts = 100;
    deltaV = velocity1 - velocity0;
    smallMove[{oldPosition_, oldVelocity_}] := (
        newPosition = oldPosition + oldVelocity*((deltaT)/(noCuts));
        newVelocity = oldVelocity + (deltaV)/(noCuts);
        {newPosition, newVelocity}
    );
    #[[1]] & /@ NestList[smallMove, {position0, velocity0}, noCuts]
);

并且重新作图和重新计算误差:

再重新计算一遍

再重新计算一遍

可以看到误差减小了,可以看到整个计算方法是没问题.一般而言,运动的时间越短,速度的改变幅度越大,那么运动的轨迹就会越陡峭,那么,区间分隔的个数就应该相应地多一些,这样呢才能把它这个误差给它控制住.

卫星的轨道半径和周期

卫星呢它也有它自己的这个速度,就假如说这个地心引力突然消失的话,这个卫星它肯定就会沿着它当时的轨道的切线方向飞出去飞到很远的地.道理呢,类似于,下雨天我们出门会打伞,如果注意到会发现,沿着伞柄旋转这个伞的时候,水滴会从从伞边飞出去溅到别人身上,就是因为我们转这个伞的时候,伞面上的水珠也相当于做着这个圆周运动,而水珠和伞面的摩擦力不足以拉着它,这个水珠它就被甩出去了,而这个使水珠在伞面上停留的力提供的就是向心力,只不过它提供的向心力达不到要求,拉不住水珠,所以水珠就飞溅出去了.

一般地,假如说一个速度方向垂直于半径的卫星的质量是$m$,速度大小是$v$,距离母星质心的距离是$r$,那么要使该卫星的运动轨道保持为一个圆,母星需要用刚好等于

$$ F_c = \frac{mv^2}{r} $$

这么大的力气来吸引这个卫星,如果说母星对卫星的引力大过这个值太多,那么卫星就会被吸到母星表面上坠毁,如果说母星对卫星的引力小过这个值,卫星的运动轨迹就会是椭圆,如果母星对卫星的引力小过这个值太多太多,那么卫星就会飞出轨道,母星就留不住这颗卫星.

一颗卫星,轨道越低,它运行速度越.为什么这么说呢?因为,我们知道,卫星维持轨道运动所需的向心力实际上是由万有引力提供的,而万有引力的大小也就取决于卫星和母星的质量,万有引力并不是因为卫星做圆周运动才有的,所以,卫星的运动状态的改变不会改变自身受到的用作向心力的万有引力的大小,也就是说:

$$ \frac{mv^2}{r} = \frac{G M m}{r^2} $$

其中的$M$是母星的质量,而$G$是万有引力常数,$r$是轨道高度,那么这个式子我们给它变一下形,就得到

$$ v = \sqrt{\frac{G M}{r}} $$

所以说,轨道高度越低,卫星的运动速度就越.我们看到低轨卫星都是十几分钟就绕地球转一圈的,而像月球这样的轨道高高的的卫星,要大约一个月才能绕着地球转一.这里的$v$表示的仅仅是卫星在单位时间内划过的弧长的长度,如果我们想知道卫星单位时间内划过的角度怎么算呢?卫星假如说要绕轨道走一圈,要走的路程有$2 \pi r$这么长,那么,卫星走一圈需要多长时间呢?那就是需要

$$ T = \frac{2 \pi r}{v} $$

这么长的时间,其中的$T$叫做周期,周期就是卫星绕母星跑一圈所需的时间,哦!这样我们同时还知道了,卫星需要长度为$T$的时间来跑完$2 \pi$的弧度,那么,单位时间内卫星划过的弧度就是

$$ \omega = \frac{2 \pi}{T} = \frac{2 \pi}{\frac{2 \pi r}{v}} = \frac{v}{r} $$

那么划过单位弧度所需的时间自然就是$\omega^{-1} = r / v$.式中的$\omega$也被叫做角速.

以上的式子在我们设计卫星时能为我们提供帮助,假如说,我们想设计一颗运行轨道为圆形的火箭,并且,我们希望卫星在长度为$T$的时间内绕地球转一圈,那么,发射的时候应该让火箭把这颗卫星推到多高的轨道上去呢?

我们可以这样想,火箭的运行速度$v$是由轨道半径$r$决定的,如果我确定了轨道半径$r$那么就等于确定了运行速度$v$,而确定了$v$和$r$就等于确定了角速度$\omega$因此也就确定了周期$T$,即速度$v$可看做是关于轨道半径$r$的一个函数

$$ v = v(r) = \sqrt{\frac{GM}{r}} $$

与此同时,周期$T$也可看做是关于$(v,r)$的函数

$$ T = T(v, r) = \frac{2 \pi r}{v} $$

也就是

$$ T = T(v, r) = T(v(r), r) = T(r) $$

也就是说周期$T$也可看做是一个仅和$r$有关的函数,并且

$$ T =\frac{2 \pi r^{3/2}}{\sqrt{G M}} $$

我们可以从这个式子中反解出$r = r(T)$,也就是反解出用周期表示的轨道高度$r$

$$ r = \frac{T^{2/3} \sqrt[3]{G M}}{\sqrt[3]{4 \pi^{2}}} $$

对于圆形轨道的卫星,我们希望卫星多久绕一圈地球,我们就可以根据上面这个式子算出要把这个卫星推到多高的位置上.

卫星运动的轨迹

我们记第$n$次迭代时,卫星的物理状态为

$$ \begin{bmatrix} {F^{(n)}} \\
{a^{(n)}} \\
{v^{(n)}} \\
{x^{(n)}} \end{bmatrix} $$

式中的$F,a,v,x$分别表示卫星受到的合力,卫星的加速度,卫星的速度和卫星的坐标,那么,卫星的物理状态的迭代实际上可表示为

$$ \begin{bmatrix} {F^{(n)}} \\
{a^{(n)}} \\
{v^{(n)}} \\
{x^{(n)}} \end{bmatrix} \rightarrow \begin{bmatrix} {F^{(n+1)}} \\
{a^{(n+1)}} \\
{v^{(n+1)}} \\
{x^{(n+1)}} \end{bmatrix} $$

更具体地,每一次迭代求新一期的物理状态的近似公式为

$$ \begin{cases} {F^{(n+1)}} \approx \frac{G M m}{d({x_c}, {x^{(n)}})^2} \frac{{x_c} - {x^{(n)}}}{|{x_c} - {x^{(n)}}|} \\
{a^{(n+1)}} = \frac{{F^{(n+1)}}}{m} \\
{v^{(n+1)}} \approx {v^{(n)}} + \Delta t_{\epsilon} {a^{(n+1)}} \\
{x^{(n+1)}} \approx {x^{(n)}} + \Delta t_{\epsilon} {v^{(n+1)}} \end{cases} $$

其中,$M$表示母星的质量,$m$表示卫星的质量,$x_c$表示母星的坐标,$x$表示卫星的坐标,$d(\cdot, \cdot)$表示欧式距离,$\Delta t_{\epsilon}$表示非常短的时间间.更加准确地,我们可以把系统第$n+1$次迭代看成是在时刻$t$完成的,而第$n$次迭代完成的世界显然就是$t-\Delta t_{\epsilon}$,那么,上式还可以写作

$$ \begin{cases} {F(t)} \approx \frac{G M m}{d({x_c}, {x(t-\Delta t_{\epsilon})})^2} \frac{{x_c} - {x(t-\Delta t_{\epsilon})}}{|{x_c} - {x(t-\Delta t_{\epsilon})}|} \\
{a(t)} = \frac{{F(t)}}{m} \\
{v(t)} \approx {v(t-\Delta t_{\epsilon})} + \Delta t_{\epsilon} {a(t)} \\
{x(t)} \approx {x(t-\Delta t_{\epsilon})} + \Delta t_{\epsilon} {v(t)} \end{cases} $$

其中$\Delta t_{\epsilon}$越小,等式左右两端的差别就越小,估计得就越准确,观察这组式子,似乎使我们想到了什么,我们先对其中的一些式子做一些等价的改写:

$$ \begin{cases} {F(t)} \approx \frac{G M m}{d({x_c}, {x(t-\Delta t_{\epsilon})})^2} \frac{{x_c} - {x(t-\Delta t_{\epsilon})}}{|{x_c} - {x(t-\Delta t_{\epsilon})}|} \\
{a(t)} = \frac{{F(t)}}{m} \\
\frac{{v(t)} - {v(t-\Delta t_{\epsilon})}}{\Delta t_{\epsilon}} \approx {a(t)} \\
\frac{{x(t)} - {x(t-\Delta t_{\epsilon})}}{\Delta t_{\epsilon}} \approx {v(t)} \end{cases} $$

观察到上式正是微分的近似形式,事实上,只要令$\Delta t_{\epsilon} \to 0$,就立刻有

$$ \begin{cases} {F(t)} = \frac{G M m}{d({x_c}, {x(t)})^2} \cdot \frac{{x_c} - {x(t)}}{|{x_c} - {x(t)}|} \\
{a(t)} = \frac{{F(t)}}{m} \\
v^{\prime}(t) = a(t) \\
x^{\prime}(t) = v(t) \end{cases} $$

对上式化简,可以得到用来描述卫星的轨迹$x(t)$的微分方程

$$ x^{\prime \prime}(t) = \frac{G M}{d({x_c}, {x(t)})^2} \cdot \frac{{x_c} - {x(t)}}{|{x_c} - {x(t)}|} $$

通过变换坐标系,我们令母星的坐标等于$0$,事实上这总是可以办到的,于是,$x(t)$的实际意义变为卫星相对母星的坐标,而上式就变为

$$ x^{\prime \prime}(t) = -\frac{G M}{|x(t)|^3} x(t) $$

求卫星的轨迹,就是从中解出$x(t)$.式中的$x$实际上是向量,我们记它的坐标形式为$x(t)=(p(t),q(t),s(t))$,那么,由以上向量形式的微分方程可以得到微分方程组

$$ \begin{cases} p^{\prime \prime} = -\frac{GM}{(p^2+q^2+s^2)^(3/2)}p \\
q^{\prime \prime} = -\frac{GM}{(p^2+q^2+s^2)^(3/2)}q \\
s^{\prime \prime} = -\frac{GM}{(p^2+q^2+s^2)^(3/2)}s \end{cases} $$

其中的$(p, q, s)$就是$t$时刻卫星相对母星的坐标.即使不会解这个微分方程组也没关系,卫星的行进轨迹还是可以估计的,用那个简单的迭代公式来估计,还可以用微分方程的数值求解方法来求解,后者相比我们自己写程序迭代求值更加推荐,因为专用的微分方程数值求解方法更加精确和高效.

gravitationalConstant = Entity[
    "PhysicalConstant", 
    "GravitationalConstant"
]["Value"] // QuantityMagnitude;

earthMass = QuantityMagnitude[Entity["Planet", "Earth"]["Mass"]];

g = gravitationalConstant;
m = earthMass;

x0 = 1000000;
y0 = 0;
z0 = 0;

v0p = 0;
v0q = 9660.48;
v0s = 0;

sol = NDSolve[
    {
        p''[t] == -(g*m*p[t])/((p[t]^2 + q[t]^2 + s[t]^2)^(3/2)),
        q''[t] == -(g*m*q[t])/((p[t]^2 + q[t]^2 + s[t]^2)^(3/2)),
        s''[t] == -(g*m*s[t])/((p[t]^2 + q[t]^2 + s[t]^2)^(3/2)),
        p'[0] == v0p,
        q'[0] == v0q,
        s'[0] == v0s,
        p[0] == x0,
        q[0] == y0,
        s[0] == z0
    },
    {
        p, q, s
    },
    {
        t, 0, 1000
    }
];

ParametricPlot3D[
    {p[t], q[t], s[t]} /. sol, 
    {t, 0, 1000}, 
    ViewPoint -> {0, 0, 2}, 
    ViewVertical -> {0, 1, 0}
]

输出结果如下图

卫星的轨迹图(垂直轨道所在平面的视角)

卫星的轨迹图(垂直轨道所在平面的视角)

微分方程真是描述连续系统演化过程的好工具!

再次回顾变速直线运动

一开始我们求物体从一个点到另一个点的运动轨迹,我们用的是迭代法来求,而现在,学习到了微分方程和差分迭代过程的联系后,我们可以把这种变速运动表示为

$$ \begin{cases} x^{\prime}(t) = v(t) \\
v^{\prime}(t) = a(t) \\
a(t) = a_0 \\
v(0) = v_0 \\
x(0) = x_0 \end{cases} $$

也就是说,物体从一个点到另一个点的变速运动的轨迹实际上也可以用微分方程来表.而卫星绕着母星的运动,其实只不过是令上列方程中的加速度$a(t)$等于母星对其造成的引力所引起的加速度$a_{c}(t)$罢了,也就是

$$ \begin{cases} x^{\prime}(t) = v(t) \\
v^{\prime}(t) = a(t) \\
a(t) = a_{c}(t) \\
v(0) = v_0 \\
x(0) = x_0 \end{cases} $$

对于万有引力下的运动,我们已经知道了

$$ \begin{cases} x^{\prime}(t) = v(t) \\
v^{\prime}(t) = a(t) \\
a(t) = -\frac{G M}{r^3} x(t) \\
v(0) = v_0 \\
x(0) = x_0 \end{cases} $$

其中$x(t)$表示卫星相对母星的坐标,而$r$表示卫星到母星的距离,多除了一次$r$是因为要把向量$x(t)$归一.而对于两点间的匀速运动,用微分方程表示起来也很简单

$$ \begin{cases} x^{\prime}(t) = v_c \\
x(0) = x_0 \end{cases} $$

请体会物理在数学中的简洁和统一!

推广和应用

下面,让我们来把我们之前写的轨道求解程序封装成一个函数,方便调用:

singleSatelliteSystem[
    masterMass_,
    gravitationalConstant_,
    relX_, 
    relY_, 
    relZ_, 
    relVX_,
    relVY_,
    relVZ_,
    timeToEvolve_
] := (

    g = gravitationalConstant;
    m = masterMass;

    x0 = relX;
    y0 = relY;
    z0 = relZ;

    v0p = relVX;
    v0q = relVY;
    v0s = relVZ;

    sol = NDSolve[
        {
            p''[t] == -(g*m*p[t])/((p[t]^2 + q[t]^2 + s[t]^2)^(3/2)),
            q''[t] == -(g*m*q[t])/((p[t]^2 + q[t]^2 + s[t]^2)^(3/2)),
            s''[t] == -(g*m*s[t])/((p[t]^2 + q[t]^2 + s[t]^2)^(3/2)),
            p'[0] == v0p,
            q'[0] == v0q,
            s'[0] == v0s,
            p[0] == x0,
            q[0] == y0,
            s[0] == z0
        },
        {
            p, q, s
        },
        {
            t, 0, timeToEvolve
        }
    ];

    sol
);

参数masterMass表示的是母星的质量,而gravitationalConstant表示的则是引力常量,relX,relY,和relZ分别表示卫星相对母星的偏移量,而relV*则表示卫星三个维度上的速度大小,timeToEvolve表示运动进行的时间长.接下来,我们顺便把单个物体匀速运动和变速运动的轨迹求解程序也写出来:

constantVelocitySystem[
    initialX_,
    initialY_,
    initialZ_,
    velocityX_,
    velocityY_,
    velocityZ_,
    timeToEvolve_
] := (

    sol = NDSolve[
        {
            x'[t] == velocityX,
            y'[t] == velocityY,
            z'[t] == velocityZ,
            x[0] == initialX,
            y[0] == initialY,
            z[0] == initialZ
        },
        { x, y, z },
        { t, 0, timeToEvolve }
    ];

    sol
);

acceleratingSystem[
    initialX_,
    initialY_,
    initialZ_,
    initialVX_,
    initialVY_,
    initialVZ_,
    accelerationX_,
    accelerationY_,
    accelerationZ_,
    timeToEvolve_
] := (

    sol = NDSolve[
        {
            x''[t] == accelerationX,
            y''[t] == accelerationY,
            z''[t] == accelerationZ,
            x'[0] == initialVX,
            y'[0] == initialVY,
            z'[0] == initialVZ,
            x[0] == initialX,
            y[0] == initialY,
            z[0] == initialZ
        },
        { x, y, z },
        { t, 0, timeToEvolve }
    ];

    sol
)

评价

本来,咱们只是想用迭代的方法,或者用近似的方法,一步步地演化出每一个离散时间点下的物体的物理状态,不曾想,哎呀!竟然就自然而然地从状态转移方程推出了微分方程来了!由此我们受到启发,状态转移方程描述的是离散系统的演化,而微分方程描述的是连续系统的演化,当时间粒度足够细时,二者的差别也足够小,这其实就是离散系统的模拟和连续系统的模拟的相似之处与底层联.

不看不知道,数学的世界真奇妙!

为了好玩physicsastronomymath

在Mathematica中实现魔方

判断一个点是否在一个多边形内部的几种方法以及它们的实现