基于GeoGebra和VPython的天體運行數(shù)值模擬*

陳海濤

(云南師范大學(xué)物理與電子信息學(xué)院 云南 昆明 650500)

金惠吉

(昆明市第八中學(xué) 云南 昆明 650222)

楊秀發(fā)

(昆明市第十四中學(xué) 云南 昆明 650106)

任 鵬

(云南師范大學(xué)實驗中學(xué) 云南 昆明 650031)

2018年4月教育部印發(fā)的《教育信息化2.0行動計劃》指出，以教育信息化支撐引領(lǐng)教育現(xiàn)代化，是新時代我國教育改革發(fā)展的戰(zhàn)略選擇，對于構(gòu)建教育強國和人力資源強國具有重要意義[1]. 2017版最新《普通高中物理課程標(biāo)準(zhǔn)》在高中物理課程的基本理念中指出，要通過多樣化的教學(xué)方式，利用現(xiàn)代信息技術(shù)，引導(dǎo)學(xué)生理解物理學(xué)的本質(zhì)，增強科學(xué)探究能力[2].

近年來，有許多學(xué)者在利用信息技術(shù)輔助物理教學(xué)方面取得了豐碩成果.例如， 2021年，文獻[3]用COMSOL軟件模擬了“磁場和磁感線”；文獻[4]用Unity3D軟件設(shè)計和仿真了光學(xué)實驗;2020年，文獻[5]將LabVIEW軟件運用到了光速測定實驗中等.

Python是一種面向?qū)ο蟮某绦蛟O(shè)計語言，語法簡潔清晰、上手容易，適合非專業(yè)計算機出身的物理教師學(xué)習(xí). VPython是建構(gòu)在Python程序設(shè)計語言之上的一個模塊，能夠方便做出三維動畫且計算功能也非常強大，很適合進行物理的仿真模擬.

GeoGebra是一款開源、免費、易上手，且支持PC端、手機端、網(wǎng)頁端等多種渠道操作的功能強大的軟件，無需編程功底就可以制作出很多物理情境的模擬.

1 牛頓的“大炮”

如圖1所示,在離地面一定高度水平拋出一物體，當(dāng)初速度小于第一宇宙速度時，物體沿橢圓曲線a、a1落地；當(dāng)速度為第一宇宙速度時，物體沿圓軌道b運行；當(dāng)初速度介于第一宇宙速度和第二宇宙速度之間時，物體沿橢圓軌道c運行；初速度等于第二宇宙速度時，物體沿拋物線軌道d離開地球，大于第二宇宙速度時物體沿雙曲線e離開地球.

圖1 牛頓的“大炮”示意圖

為了模擬這一現(xiàn)象，可先由牛頓第二定律和幾何關(guān)系，建立微分方程組

1.1 Vpython數(shù)值模擬方法

將上述微分方程和初始條件轉(zhuǎn)化為計算機代碼，基于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)

#參數(shù)設(shè)置

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=[] #衛(wèi)星

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#計算衛(wèi)星瞬

時徑向矢量

Fn=G*m_earth*m_satellite/(mag(axis))**2#計算萬有引力數(shù)值大小

FG=Fn*axis.norm() #計算萬有

引力矢量

b.a=FG/m_satellite #計算衛(wèi)星加

速度矢量

b.v+=b.a*dt#計算衛(wèi)星速度矢量

b.pos += b.v*dt#計算衛(wèi)星位置矢量

運行效果如圖2所示．

圖2 Python運行效果圖

接著還可以使用Vpython中的繪圖語句實時繪制各衛(wèi)星在軌道上的能量變化情況，如圖3～5所示.

圖3 衛(wèi)星在橢圓軌道運行時的能量變化情況

圖 4 衛(wèi)星在拋物線軌道上運行時的能量變化情況

圖 5 衛(wèi)星在雙曲線軌道上運行時的能量變化情況

1.2 GeoGebra數(shù)值模擬方法

在GeoGebra中，無需編程，基于解常微分方程組的指令，僅需要非常簡單的步驟就可以模擬天體運動.

首先在GeoGebra代數(shù)區(qū)中定義參數(shù)滑動條，為方便起見，可設(shè)G=1，M=10，m=1，環(huán)繞天體初始位置的橫坐標(biāo)x_{01}=10，縱坐標(biāo)y_{01}=0，初速度水平分量v_{x01}=0，豎直分量v_{y01}=1，接著根據(jù)微分方程分別輸入

再根據(jù)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}”.指令中第一部分是導(dǎo)數(shù)列表，第二部分的“0”為t的初始值，第三部分分別是t=0時的“x_{1}”“y_{1}”“v_{x1}”和“v_{y1}”的值，第四部分的1000為t的末值.再輸入“l(fā)en=length(numericalIntegral1)”，給出構(gòu)成圖像的點的數(shù)目；輸入“t′=Slider(0,1,1/len)”， 其中“1/len”為增量，這樣t′每次增加都會對應(yīng)于numericalIntegral1中的下一個點.接著設(shè)中心天體所在位置“A=(0,0)”，為了得到環(huán)繞天體的位置，先輸入“B=Point(numericalIntegral1,t′)”，得到“(t,x_{1})”點，該描點指令中第二個t′是路徑值，取值在0～1之間，決定了取值在numerical-Integral1中的相應(yīng)位置，例如若取“0”則為第一個點，“1”則為最后一個點，“1/len”則為第二個點.再輸入“C=Point(numericalIntegral1_{2},t′)”得到“(t,y_{1})”點，輸入“D=(y(B),y(C))”即可得到環(huán)繞天體位置坐標(biāo).最后，啟動t′動畫即可觀察到環(huán)繞天體的運動，t′乘上1000即為真實時間t.

改變不同的環(huán)繞天體初速度大小即可觀察到如圖6～8所示的軌跡圖，可以發(fā)現(xiàn)當(dāng)初速度正好等于臨界速度1時，軌跡的確是一個圓且環(huán)繞一周用時和理論計算是一致的，當(dāng)初速度大小為1.3時，可以觀察到軌跡的確是一個橢圓，且環(huán)繞一周大約用了364 s，當(dāng)初速度大小為2.7 m/s，即大于第二宇宙速度時，可以觀察到環(huán)繞天體脫離地球的束縛.

圖6 設(shè)置“v_{y01}=1”時

圖7 設(shè)置“v_{y01}=1.3”時

圖8 設(shè)置“v_{y01}=2.7”時

2 拉格朗日點

為了展示出“太陽-地球”雙星系統(tǒng)的運行情況，可設(shè)G=6.67，太陽質(zhì)量ms為100，地球質(zhì)量me為30，衛(wèi)星質(zhì)量為m，太陽到質(zhì)心C的距離為6，質(zhì)心C位于原點.根據(jù)質(zhì)心公式即可計算出地球到質(zhì)心C的距離為20.根據(jù)理論易推導(dǎo)出L4和L5與太陽和地球的位置嚴(yán)格構(gòu)成等邊三角形的關(guān)系[6～8]，如圖9所示.

圖9 L4和L5拉格朗日點示意圖

根據(jù)牛頓第二定律和對稱性即可得出各天體的微分方程

2.1 Vpython數(shù)值模擬方法

將拉格朗日點的數(shù)學(xué)模型和初始參數(shù)轉(zhuǎn)化為計算機代碼，基于Vpython庫用Python編程如下(部分代碼)：

t=0;dt=0.001;m_earth=30;m_sun=100;G=6.67 #設(shè)置參數(shù)

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)#質(zhì)心

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)

運行結(jié)果如圖10所示.

圖10 L4和L5拉格朗日點Python模擬效果圖

2.2 GeoGebra數(shù)值模擬方法

在GeoGebra中，類似地，基于解常微分方程組的指令，無需編程就可模擬出兩衛(wèi)星分別處在拉格朗日點L4和L5時太陽、地球以及兩衛(wèi)星的運動.

首先在GeoGebra代數(shù)區(qū)中定義參數(shù)滑動條，為方便起見，可設(shè)G=6.67，太陽質(zhì)量m1=100，地球質(zhì)量m2=30，兩衛(wèi)星質(zhì)量由于遠遠小于地球和太陽質(zhì)量，不妨分別設(shè)為m3=0，m4=0，接著，若讓質(zhì)心處在原點(0,0)處，則根據(jù)理論推導(dǎo)易得太陽、地球、處在L4位置和L5位置的衛(wèi)星的橫坐標(biāo)應(yīng)分別定義為:

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

縱坐標(biāo)分別定義為：

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})))

接著根據(jù)微分方程輸入

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)=

根據(jù)對稱性很容易給出另外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 結(jié)束語

本文用Vpython和GeoGebra分別成功模擬了牛頓的“大炮”以及拉格朗日點的運行軌跡，其中前者采用了歐拉法的數(shù)值模擬方法，后者則采用了GeoGebra模擬物理情境方面的最新技術(shù)——解常微分方程組指令，兩者都通過物理情境模擬生動展示出了理論的預(yù)測結(jié)果，不僅驗證了理論，同時也加深了我們對理論的理解.若將這個技術(shù)運用到物理的教學(xué)與學(xué)習(xí)中，能夠使自身以及學(xué)生更好地理解數(shù)學(xué)模型、計算機算法和相應(yīng)物理現(xiàn)象之間的聯(lián)系，提高跨學(xué)科的綜合運用能力．