陳海濤

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

金惠吉

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

楊秀發(fā)

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

任 鵬

(云南師范大學(xué)實(shí)驗(yàn)中學(xué) 云南 昆明 650031)

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

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

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

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

1 牛頓的“大炮”

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

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

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

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

將上述微分方程和初始條件轉(zhuǎn)化為計(jì)算機(jī)代碼，基于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#計(jì)算衛(wèi)星瞬

時(shí)徑向矢量

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

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

引力矢量

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

速度矢量

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

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

運(yùn)行效果如圖2所示．

圖2 Python運(yùn)行效果圖

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

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

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

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

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

在GeoGebra中，無需編程，基于解常微分方程組的指令，僅需要非常簡(jiǎn)單的步驟就可以模擬天體運(yùn)動(dòng).

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

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

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

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

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

2 拉格朗日點(diǎn)

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

圖9 L4和L5拉格朗日點(diǎn)示意圖

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

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

將拉格朗日點(diǎn)的數(shù)學(xué)模型和初始參數(shù)轉(zhuǎn)化為計(jì)算機(jī)代碼，基于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)

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

圖10 L4和L5拉格朗日點(diǎn)Python模擬效果圖

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

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

首先在GeoGebra代數(shù)區(qū)中定義參數(shù)滑動(dòng)條，為方便起見，可設(shè)G=6.67，太陽質(zhì)量m1=100，地球質(zhì)量m2=30，兩衛(wèi)星質(zhì)量由于遠(yuǎn)遠(yuǎn)小于地球和太陽質(zhì)量，不妨分別設(shè)為m3=0，m4=0，接著，若讓質(zhì)心處在原點(diǎn)(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ù)對(duì)稱性很容易給出另外3個(gè)天體的微分方程，最后得到有16個(gè)微分方程的微分方程組.接著再輸入解微分方程組指令解出16個(gè)微分方程：

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中例子類似，不再贅述.最終運(yùn)行效果如圖11所示.

圖11 L4和L5拉格朗日點(diǎn)GeoGebra模擬效果圖

3 結(jié)束語

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