基于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模擬物理情境方面的最新技術——解常微分方程組指令，兩者都通過物理情境模擬生動展示出了理論的預測結果，不僅驗證了理論，同時也加深了我們對理論的理解.若將這個技術運用到物理的教學與學習中，能夠使自身以及學生更好地理解數學模型、計算機算法和相應物理現象之間的聯系，提高跨學科的綜合運用能力．