基于GeoGebra和VPython的天体运行数值模拟*

陈海涛

(云南师范大学物理与电子信息学院 云南 昆明 650500)

金惠吉

(昆明市第八中学 云南 昆明 650222)

杨秀发

(昆明市第十四中学 云南 昆明 650106)

任 鹏

(云南师范大学实验中学 云南 昆明 650031)

2018年4月教育部印发的《教育信息化2.0行动计划》指出，以教育信息化支撑引领教育现代化，是新时代我国教育改革发展的战略选择，对于构建教育强国和人力资源强国具有重要意义[1]. 2017版最新《普通高中物理课程标准》在高中物理课程的基本理念中指出，要通过多样化的教学方式，利用现代信息技术，引导学生理解物理学的本质，增强科学探究能力[2].

近年来，有许多学者在利用信息技术辅助物理教学方面取得了丰硕成果.例如， 2021年，文献[3]用COMSOL软件模拟了“磁场和磁感线”；文献[4]用Unity3D软件设计和仿真了光学实验;2020年，文献[5]将LabVIEW软件运用到了光速测定实验中等.

Python是一种面向对象的程序设计语言，语法简洁清晰、上手容易，适合非专业计算机出身的物理教师学习. VPython是建构在Python程序设计语言之上的一个模块，能够方便做出三维动画且计算功能也非常强大，很适合进行物理的仿真模拟.

GeoGebra是一款开源、免费、易上手，且支持PC端、手机端、网页端等多种渠道操作的功能强大的软件，无需编程功底就可以制作出很多物理情境的模拟.

1 牛顿的“大炮”

如图1所示,在离地面一定高度水平抛出一物体，当初速度小于第一宇宙速度时，物体沿椭圆曲线a、a1落地；当速度为第一宇宙速度时，物体沿圆轨道b运行；当初速度介于第一宇宙速度和第二宇宙速度之间时，物体沿椭圆轨道c运行；初速度等于第二宇宙速度时，物体沿抛物线轨道d离开地球，大于第二宇宙速度时物体沿双曲线e离开地球.

图1 牛顿的“大炮”示意图

为了模拟这一现象，可先由牛顿第二定律和几何关系，建立微分方程组

1.1 Vpython数值模拟方法

将上述微分方程和初始条件转化为计算机代码，基于Vpython库用Python编程如下：

from vpython import*

scene=canvas(width=1200,height=1000,background=vector(1,1,1))

Re=6.4*10**6;H=1.5*Re;m_earth=6*10**24;m_satellite=1000;G=6.67*10**

(-11)

#参数设置

v0=(G*m_earth/H)**0.5

T=2*pi*H/v0

earth=sphere(pos=vector(0,0,0),radius=Re,texture=textures.earth,opacity=0.6) #地球

satellite=[] #卫星

for i in range(5,18,1):

satellite.append(sphere(pos=vector(0,H,0),v=vector(0.1*i*v0,0,0),radius=

0.09*Re,a=vector(0,0,0),color=color.red,make_trail=1) )

t=0;dt=0.01

while t<=T:

t =t+ dt

rate(10000)

for b in satellite:

axis=earth.pos - b.pos#计算卫星瞬

时径向矢量

Fn=G*m_earth*m_satellite/(mag(axis))**2#计算万有引力数值大小

FG=Fn*axis.norm() #计算万有

引力矢量

b.a=FG/m_satellite #计算卫星加

速度矢量

b.v+=b.a*dt#计算卫星速度矢量

b.pos += b.v*dt#计算卫星位置矢量

运行效果如图2所示．

图2 Python运行效果图

接着还可以使用Vpython中的绘图语句实时绘制各卫星在轨道上的能量变化情况，如图3～5所示.

图3 卫星在椭圆轨道运行时的能量变化情况

图 4 卫星在抛物线轨道上运行时的能量变化情况

图 5 卫星在双曲线轨道上运行时的能量变化情况

1.2 GeoGebra数值模拟方法

在GeoGebra中，无需编程，基于解常微分方程组的指令，仅需要非常简单的步骤就可以模拟天体运动.

首先在GeoGebra代数区中定义参数滑动条，为方便起见，可设G=1，M=10，m=1，环绕天体初始位置的横坐标x_{01}=10，纵坐标y_{01}=0，初速度水平分量v_{x01}=0，竖直分量v_{y01}=1，接着根据微分方程分别输入

再根据GeoGebra的解常微分方程组指令输入“NSolveODE({x_{1}′,y_{1}′,v_{x1}′,v_{y1}′},0,{x_{01},y_{01},v_{x01},v_{y01}},1000)”，直接输入后，软件会自动分别给出“x_{1}-t”“y_{1}-t”“v_{x1}-t”“v_{y1}-t”图像，自动命名为“numericalIntegral1、numericalIntegral1_{2}、numericalIntegral1_{3} 、numericalIntegral1_{4}”.指令中第一部分是导数列表，第二部分的“0”为t的初始值，第三部分分别是t=0时的“x_{1}”“y_{1}”“v_{x1}”和“v_{y1}”的值，第四部分的1000为t的末值.再输入“len=length(numericalIntegral1)”，给出构成图像的点的数目；输入“t′=Slider(0,1,1/len)”， 其中“1/len”为增量，这样t′每次增加都会对应于numericalIntegral1中的下一个点.接着设中心天体所在位置“A=(0,0)”，为了得到环绕天体的位置，先输入“B=Point(numericalIntegral1,t′)”，得到“(t,x_{1})”点，该描点指令中第二个t′是路径值，取值在0～1之间，决定了取值在numerical-Integral1中的相应位置，例如若取“0”则为第一个点，“1”则为最后一个点，“1/len”则为第二个点.再输入“C=Point(numericalIntegral1_{2},t′)”得到“(t,y_{1})”点，输入“D=(y(B),y(C))”即可得到环绕天体位置坐标.最后，启动t′动画即可观察到环绕天体的运动，t′乘上1000即为真实时间t.

改变不同的环绕天体初速度大小即可观察到如图6～8所示的轨迹图，可以发现当初速度正好等于临界速度1时，轨迹的确是一个圆且环绕一周用时和理论计算是一致的，当初速度大小为1.3时，可以观察到轨迹的确是一个椭圆，且环绕一周大约用了364 s，当初速度大小为2.7 m/s，即大于第二宇宙速度时，可以观察到环绕天体脱离地球的束缚.

图6 设置“v_{y01}=1”时

图7 设置“v_{y01}=1.3”时

图8 设置“v_{y01}=2.7”时

2 拉格朗日点

为了展示出“太阳-地球”双星系统的运行情况，可设G=6.67，太阳质量ms为100，地球质量me为30，卫星质量为m，太阳到质心C的距离为6，质心C位于原点.根据质心公式即可计算出地球到质心C的距离为20.根据理论易推导出L4和L5与太阳和地球的位置严格构成等边三角形的关系[6～8]，如图9所示.

图9 L4和L5拉格朗日点示意图

根据牛顿第二定律和对称性即可得出各天体的微分方程

2.1 Vpython数值模拟方法

将拉格朗日点的数学模型和初始参数转化为计算机代码，基于Vpython库用Python编程如下(部分代码)：

t=0;dt=0.001;m_earth=30;m_sun=100;G=6.67 #设置参数

earth=sphere(pos=vector(20,0,0),radius=1,texture=textures.earth,make_trail=1)#地球

sun=sphere(pos=vector(-6,0,0),radius=2,color=color.red,make_trail=1) #太阳

x0=sphere(pos=vector(0,0,0),radius=0.1,color=color.red)#质心

r=26;r1=20;r2=6;Fn=G*m_earth*m_

sun/(r**2);w=sqrt((G*(m_earth+m_sun))/

(r**3));period=2*pi/w

earth_v0=(G*m_sun*r1/(r**2))**0.5;earth.v=vector(0,earth_v0,0)

sun_v0=(G*m_earth*r2/(r**2))**0.5;sun.v=vector(0,-sun_v0,0)

m_prober4=1;prober4=sphere(pos=vector(r*cos(pi/3)-6,r*sin(pi/3),0),radius=0.5,color=color.red,make_trail=1)

r4=sqrt(r2**2+r**2-2*r2*r*cos(pi/3));cos_beta=(26**2+r4**2-r2**2)/(2*26*r4)

β1=acos(cos_beta);α1=(pi/6-β1);v4_0=w*(r4)

v4=vector(-v4_0*cos(α1),v4_0*sin(α1),0)

m_prober5=1;prober5=sphere(pos=vector(r*cos(pi/3)-6,-r*sin(pi/3),0),radius=0.5,color=color.red,make_trail=1)

r5=sqrt(r2**2+r**2-2*r2*r*cos(pi/3));cos_beta=(26**2+r5**2-r2**2)/(2*26*r5)

β2= acos(cos_beta); α2=(pi/6-β2);v5_0=w*(r5)

v5=vector(v5_0 *cos(α2),v5_0 *sin(α2),0)

运行结果如图10所示.

图10 L4和L5拉格朗日点Python模拟效果图

2.2 GeoGebra数值模拟方法

在GeoGebra中，类似地，基于解常微分方程组的指令，无需编程就可模拟出两卫星分别处在拉格朗日点L4和L5时太阳、地球以及两卫星的运动.

首先在GeoGebra代数区中定义参数滑动条，为方便起见，可设G=6.67，太阳质量m1=100，地球质量m2=30，两卫星质量由于远远小于地球和太阳质量，不妨分别设为m3=0，m4=0，接着，若让质心处在原点(0,0)处，则根据理论推导易得太阳、地球、处在L4位置和L5位置的卫星的横坐标应分别定义为:

x_{10}=-6、x_{20}=20、x_{30}=7、x_{40}=7

纵坐标分别定义为：

y_{10}=0、y_{20}=0、y_{30}= 26 sin(π/3)、y_{40}=-26 sin(π/3)

初速度水平分量分别定义为：

v_{x10}=0、v_{x20}=0、v_{x30}=-ω r_{4}

cos(α)、v_{x40}= ω r_{4} cos(α)

初速度竖直分量分别定义为：

v_{y10}=-sqrt(G*30 ((6)/(26^(2))))、v_{y20}=sqrt(G*100 ((20)/(26^(2))))、v_{y30}=ω r_{4} sin(α)、v_{y40}=ω r_{4} sin(α).其中，ω= sqrt(G ((m1+m2)/(( x_{20}-x_{10})^3)))，r_{4}=sqrt(6^(2)+26^(2)-2*26*6 ((1)/(2))),α=((π)/(6))-cos^(-1)(((26^(2)+r_{4}^(2)-6^(2))/(2*26 r_{4})))

接着根据微分方程输入

vy1,vx2,vy2,vx3,vy3,vx4,vy4)=vx1

vy1,vx2,vy2,vx3,vy3,vx4,vy4)=vy1

vy1,vx2,vy2,vx3,vy3,vx4,vy4)=

vy1,vx2,vy2,vx3,vy3,vx4,vy4)=

根据对称性很容易给出另外3个天体的微分方程，最后得到有16个微分方程的微分方程组.接着再输入解微分方程组指令解出16个微分方程：

NSolveODE({x_{1}′,y_{1}′,x_{2}′,y_{2}′,x_{3}′,y_{3}′,x_{4}′,y_{4}′,vx_{1}′,vy_{1}′,vx_{2}′,vy_{2}′,vx_{3}′,vy_{3}′,vx_{4}′,vy_{4}′},0,{x_{10},y_{10},x_{20},y_{20},x_{30},y_{30},x_{40},y_{40},v_{x10},v_{y10},v_{x20},v_{y20},v_{x30},v_{y30},v_{x40},v_{y40}},100)

再接着操作步骤和1.2中例子类似，不再赘述.最终运行效果如图11所示.

图11 L4和L5拉格朗日点GeoGebra模拟效果图

3 结束语

本文用Vpython和GeoGebra分别成功模拟了牛顿的“大炮”以及拉格朗日点的运行轨迹，其中前者采用了欧拉法的数值模拟方法，后者则采用了GeoGebra模拟物理情境方面的最新技术——解常微分方程组指令，两者都通过物理情境模拟生动展示出了理论的预测结果，不仅验证了理论，同时也加深了我们对理论的理解.若将这个技术运用到物理的教学与学习中，能够使自身以及学生更好地理解数学模型、计算机算法和相应物理现象之间的联系，提高跨学科的综合运用能力．