求解物理方程
——Matlab軟件在物理教學中的應用之一

朱國強

(紹興市第一中學 浙江 紹興 312000)

在1859年，中國清代數學家李善蘭開始將英文“equation”——指含有未知數的等式——一詞翻譯為“方程”．

方程按照學術程度可以分為初等數學方程和高等數學方程．方程按照解的形式可以分為數值方程和符號方程．方程按照學科分類可以分為代數方程、幾何方程、三角方程、反三角方程、邏輯方程、微分方程、積分方程等．例如中學生熟悉的初等數學領域的代數學中的方程，可稱為初等代數方程，包括整式、分式方程，一元、多元方程，有理式、無理式方程等．

兩個或兩個以上的方程的組合叫做方程組．解方程組包括數值方程組運算和符號方程組運算．Matlab不僅僅以數值運算見長，而且使用了Maple的內核，符號運算也極其厲害．見識過Matlab解符號方程組后，許多中學生深感驚奇．

1 符號變量和表達式

(1)定義符號變量

Matlab提供了函數syms，一次可以定義多個符號變量．其一般格式為：syms arg1 arg2 …argn．

(2)符號表達式

2 Solve函數及用法

Matlab內置Solve函數提供了初等數學方程的求解方案．Solve函數隸屬于符號數學工具箱(Symbolic Math Toolbox)．其用法為：

g = solve(eq)

g = solve(eq,var)

g = solve(eq1,eq2，..., ，eqn)

g = solve(eq1，eq2，...，eqn，var1，var2，...，varn)

其中，eqn為第n個方程的符號表達式，varn為第n個變量．

2.1 求解數值方程

應用Solve函數解方程通常給出的是解析解．想要獲得方程的數值解，有多種方法．第一種是應用subs函數，subs函數用于替換求值．第二種是應用eval函數，計算一個表達式的值．分別舉例說明．

(1)Solve求出解析解，subs求出數值解

Subs函數的常見格式為R = subs(S) ，用調用函數的值或者Matlab工作空間中的值來替代符號表達式S中出現的所有變量．

R = subs(S， new)，用new來代替表達式S中的默認變量；

R = subs(S，old，new) ，用new符號變量或數值或表達式來替代來符號表達式S中old符號變量或字符串表達式．

【例1】質量為m=1 kg的金屬塊放在水平桌面上，在與水平方向斜向上成θ=37°，大小為F=10 N的拉力作用下，向右做直線運動．已知金屬塊與桌面間的動摩擦因數為μ=0.5．重力加速度g=10 m/s2．求金屬塊的加速度大小．

圖1

解析：該題屬于應用牛頓運動定律解決問題中的已知受力情況求運動情況的類型．主要方程為

Fcosθ-f=ma

Fsinθ+N-mg=0

f=μN

第一步，建立模型，列方程組，求出解析解．

輸入?syms F f N miu thita m g a;

S=solve(′F*cos(thita)-f=m*a′，′F*sin(thita)+N=m*g′，′f=miu*N′，′a′，′f′，′N′);

a=simple(S.a)

解得a =(F*cos(thita)+miu*F*sin(thita)-miu*m*g)/m，其中simple函數用來簡化表達式．

即Matlab解得

第二步，根據題意，代入數據，得到數值解．

F=10;miu=0.5;thita=37*pi/180;g=10;

m=1;

a=subs(a)

解得a=5.995 4，即加速度大小為5.995 4 m/s2．

(2)Solve求出解析解，eval求出數值解

【例2】(2013年高考江蘇卷第13題)

如圖2所示，勻強磁場中有一矩形閉合線圈abcd，線圈平面與磁場垂直． 已知線圈的匝數N=100，邊長ab=1.0 m，bc=0.5 m，電阻r=2 Ω． 磁感應強度B在0～1 s內從零均勻變化到0.2 T．在1～5 s內從0.2 T均勻變化到-0.2 T，取垂直紙面向里為磁場的正方向．求：

(ⅰ)0.5 s時線圈內感應電動勢的大小E；

(ⅱ)在1～5 s內通過線圈的電荷量q；

(ⅲ)在0～5 s內線圈產生的焦耳熱Q．

圖2

解析：

該題知識點為法拉第電磁感應定律

電荷量

焦耳熱

Q=I2Rt

第一步，根據題意，建立變量，并賦數值．

syms L1 L2 r N t1 t2 t3 B0 B1 B2 E q Q

L1=1;L2=0.5;r=2;N=100;B0=0;B1=

0.2;B2=-0.2;t1=0.5;t2=1;t3=5;

第二步，建立模型，列方程組，求解方程．

A=solve(′E=N*(B1-B0)*L1*L2/(t2-0)′，′q=N*(B1-B2)*L1*L2/r′，′Q=(N*(B1-

E=eval(A.E)

q=eval(A.q)

Q=eval(A.Q)

解得E=10；q=10；Q=100.000 0．即感應電動勢的大小E=10 V ，電荷量q=10 C ，焦耳熱Q=100 J．

2.2 求解符號方程

基本方法是solve(eq1，eq2，…，eqn，v1，v2，…，vn)，即求解由表達式eq1，eq2，…，eqn組成的方程組，求解變量分別為v1，v2，…，vn．

【例3】(2008年高考全國Ⅰ卷第23題)

已知O，A，B，C為同一直線上的4點，AB間的距離為L1，BC間的距離為L2，一物體自O點靜止起出發，沿此直線做勻加速運動，依次經過A，B，C三點．已知物體通過AB段與通過BC段所用時間相等，求O與A的距離．

解法一：(人工求解)

設物體加速度為a，經過OA段時間為t，經過AB段和BC段的時間均為T，O與A間的距離為x．

對物體，OA段

(1)

(2)

(3)

聯立式(1)～(3)，解得

解法二：(計算機求解)

輸入syms x a t T L1 L2 v1;

A.x

本題屬于方程個數少于未知數個數的情形，有5種以上的解法．教學實踐表明，一些學生列出方程組、求解，需要相當長的時間，有一定的難度，而用Matlab求解則無此困難．

通過以上例子可以看出，用Matlab解各類初等數學方程組是非常方便的，只要注意正確書寫參數就可以了．

3 求解高等數學方程

高中物理教師在教研或者指導學生競賽，需要用到高等數學來分析解決問題．一般常用的有積分、微分運算，求解微分方程等．Matlab對問題的描述和求解符合人們的思維方式和數學表達習慣，可以快速地求解高等數學方程．

3.1 積分函數int

函數int的用法為S=int (y，v，a，b)，其中y是被積函數，v是積分變量，a與b積分下、上限．不定積分的輸出結果是符號表達式，定積分的輸出結果是符號表達式或數值．

【例4】計算純電阻電路中正弦交流電i=Imsinωt在一個周期內的平均功率和電流的有效值．電阻R已知．

解析：

瞬時功率p=i2R=Im2Rsin2ωt，一個周期內的平均功率

運用Matlab求解的程序如下：

syms Im omiga t R I

i=Im*sin(omiga*t);

(2*pi/omiga)

解得

即

又將p代入

I=simple(I)

解得

即有效值

3.2 dsolve函數解常微分方程

dsolve函數用法：S = dsolve(eqn,cond)，用初始或邊界條件cond求解常微分方程eqn．

【例5】如圖3所示，平行導軌上放置一金屬桿AB，質量為m，長為L，在導軌的上端接有電阻R．勻強磁場B垂直導軌平面向里．當AB桿以初速度v0向運動時，求：

(1)AB桿能夠移動的距離；

(2)在移動過程中電阻R上放出的焦耳熱．

圖3

解析：

輸入?syms v t B L R m v0

解得

即方程的特解為

輸入?syms v t B L R m v0 x

x=simple(x)

解得

即方程的特解為

當時間t→∞時，桿移動距離為

(2)桿在移動過程中產生的焦耳熱元為

整個運動過程中產生的焦耳熱為

輸入syms x v a B L m t v0 R

解得

【例6】如圖4所示，一個矩形的金屬線框，邊長分別為a和b(b足夠長)．金屬線框的質量為m，自感系數為L，忽略電阻．線框的長邊與x軸平行，它以速度v0沿x軸的方向從磁場外進入磁感應強度為B0的均勻磁場中，B0的方向垂直矩形線框平面．求矩形線框在磁場中速度與時間的關系式v=v(t)和沿x軸方向移動的距離與時間的關系式x=x(t)．

圖4

解析：

由于b邊很長，所以線框只有右邊框切割磁感線，當線框速度為v時，產生的動生電動勢ε=B0av；當線框中的電流為i時，產生的自感電動勢的大小為

根據歐姆定律得

ε+εL=iR

由于不計電阻，所以有

(4)

右邊框受力F=iaB0，根據牛頓第二定律，得

求導得

(5)

聯立式(4)、(5)得微分方程

初始條件當t=0時，v=v0，a=0．

輸入>>syms x v a B0 L m t v0

解得

即速度方程為

又由于

所以

初始條件t=0時，x=0．

輸入syms x v a B0 L m t v0

(1/2)*t)，t)

解得

即位移方程為

參考文獻

1 黃忠霖，黃京．Matlab 符號運算及其應用．北京：國防工業出版社，2004

2 朱國強．應用Matlab求解物理方程．中國教育信息化，2012(3)

3 朱國強．高中物理培優基礎教程(第二版)．杭州：浙江大學出版社，2012,318～319，441～442