Mathematica在分離變量法中的應用及極限情況下的討論*

汪夢雅 曹 蓉 夏杰楨 吳 琪

(西藏大學理學院 西藏 拉薩 850000)

靜電場問題是電磁學理論之后又一個重要的理論知識，是電動力學里的一個核心問題[1,2].其基本問題就是求解關于靜電勢的泊松方程的邊值問題，解法主要有分離變量法、格林函數法、鏡像法和電多極矩法等，其中最常用的兩種解法是分離變量法和鏡像法[3].

由于電動力學的數學要求遠比普通理工科專業高，所以為了能將學生從復雜的微分方程、偏微分方程的求解過程中解放出來，引入了Mathematica[4].Mathematica是美國WolframResearch公司開發的數學軟件，是現在使用最廣的軟件之一.該軟件除了很好地結合了數值和符號計算之外，還具有強大的繪圖功能，操作簡單，界面友好，代碼形式接近于自然語言，廣泛用于數學、物理、生物等領域[5～7].

本文利用Mathematica求解靜電場中分離變量法的問題出發，并討論在這種問題下的一些極限情況來探討Mathematica在靜電場中的應用.通過Mathematica的使用消除學生對復雜公式的畏懼，使學生可以擺脫繁瑣的數學計算，并且通過Mathematica的繪圖使學生對物理圖像有個較為直觀準確的認識，激發學生的興趣，促進課程的學習.

1 靜電場的標勢及其微分方程

1.1 靜電場的標勢

眾所周知，相距為dl的兩點的電勢差為[8]

dφ=-E·dl

(1)

然后由于

(2)

因此，電場強度E等于電勢φ的負梯度

(3)

1.2 靜電場的微分方程和邊值關系

研究一個電荷對與它鄰近的電場是怎樣作用的，一點上的電場和它鄰近的電場又是怎么聯系的，即要找出電荷和電場相互作用規律的微分形式，而在導體表面或其他邊界上場和電荷的相互關系則由邊界條件反映出來.這種問題稱為邊值問題，即求微分方程滿足邊界條件的解.

在均勻同向的介質中，D=εE，再把式(3)代入

(4)

其中ρ是自由電荷密度.式(4)是靜電勢滿足的基本微分方程，也就是泊松方程，給出邊界條件就可以確定電勢φ的解.

在兩介質界面上，靜電勢所滿足的邊值關系分別為

φ1=φ2

(5)

(6)

以上是邊值關系的一般形式，但是在靜電問題中，常常有一些導體存在，由于導體的特殊性質，在導體表面上的邊值關系有它的特點.設導體表面所帶自由電荷面密度為σ，它外面的介質電容率為ε，則由式(5)、(6)和導體靜電條件得到導體表面的邊界條件

φ=常量

(7)

(8)

靜電學的基本問題就是求出在每個均勻區域內滿足泊松方程，在所有分界面上滿足邊值關系和在所研究的整個區域邊界上滿足邊界條件的電勢的解.

2 分離變量法求解的方程和通解

在許多實際問題中，靜電場是由帶電導體決定的[9].比如，電容器內部的電場是由作為電極的兩個導體板上所帶電荷決定的，電子光學系統的靜電透鏡內部，電場是由分布于電極上的自由電荷決定的.這些問題都有一個共同的特點，就是自由電荷只出現在一些導體的表面上，而在空間中沒有其他自由電荷分布.這時，當選擇這些導體表面作為區域V的邊界，則在V內部自由電荷密度ρ=0.此時泊松方程可以化為較為簡單的拉普拉斯方程(拉式方程)，即

(9)

因此，這一類問題的解法就是求滿足邊界條件的拉普拉斯方程的解.此拉普拉斯方程的通解可以用分離變量法求出，在求解前我們要先選擇合適的坐標系，常用的坐標系有球坐標系和柱坐標系，在這里我們寫出用球坐標系得出的通解形式，球坐標用(r,θ,φ)表示，其中r為點到原點的距離，θ為極角，φ為方位角，拉式方程在球坐標中的通解為

(10)

(11)

其中Pn(cosθ)為勒讓德函數，an和bn是任意常數，由邊界條件或者說由邊值關系確定.

在每一個沒有電荷分布的區域內，φ滿足拉普拉斯方程，其通解已由式(10)和式(11)給出，剩下的問題就是由邊界條件確定這些通解中所含的任意常數，得到滿足邊界條件的特解.

3 分離變量法的實際應用

【例題】一半徑為a，介電常數為ε的介質球放置在均勻電場E0中，如圖1所示，求介質球內、外的電位及電場.

圖1 題圖

分析：介質球在外電場中極化，在它表面上產生束縛電荷.這些束縛電荷激發的電場疊加在原外電場E0上，得總電場E.束縛電荷分布和總電場E互相制約，邊界條件正確的反映這種制約關系.

球的半徑為a，球外為真空，這問題具有軸對稱性，對稱軸為通過球心沿外電場E0方向的軸線，取此軸為極軸.介質球的存在使空間分為兩均勻區域，球外區域和球內區域，兩區域內部沒有自由電荷，因此電勢φ都滿足拉普拉斯方程，以φ1代表球外區域的電勢，φ2代表球內的電勢，由式(11)得出的兩區域的通解為

(12)

(13)

式中an,bn,cn,dn是待定常數.

邊界條件有：

(1)無窮遠處，E→E0，均勻電場中空間任一點P的電勢φ(P)=φ0-E0·x，由于均勻電場可以看作由無窮大平行板電容器產生，其電荷分布不在有限區域內，因此不能選φ(∞)=0.若選φ0=0，則有φ(P)=-E0·x，可得

φ1→-E0rcosθ=-E0rP1(cosθ)

(14)

所以有

a1=-E0an=0(n≠1)

(15)

(2)r=0處，φ2應為有限值，因此

dn=0

(16)

(3)在介質球面上(r=a)

(17)

把式(12)和(13)代入得

(18)

比較P1的系數得

(19)

由式(19)解出

(20)

比較式(18)其他Pn項的系數可解出

bn=cn=0n≠1

(21)

所有常數已經定出，因此本問題的解為

(22)

(23)

為了簡化運算，用Mathematica來求解本題，并畫出介質球附近的電場線和等位面圖.所需代碼如下所示：

Clear["Global′*"]

v1=Sum[(An*R^n+Bn/R^(n+1))*Pn,{n,0,m})]

v2=Sum[(Cn*R^n+Dn/R^(n+1))*Pn,{n,0,m})]

v1=Sum[(An*R^n+Bn/R^(n+1))*Pn,{n,0,1})]/.{P0→1,P1→Cos[θ]}

v2=Sum[(Cn*R^n+Dn/R^(n+1))*Pn,{n,0,1})]/.{P0→1,P1→Cos[θ]}

(v1/.{B0→0,B0→0})==-E0*R*Cos[θ]

A0=0;A1=-E0;

D0=0;D1=0;

(v1/.R→a)==(v2/.R→a)

(ε0D[v1,R]/.R→a)==(εD[v2,R]/.R→a)

B0=0;

C0=0;

sol=Solve[{eq1,eq3},{B1,C1}]

v11=v1/.sol//First

v22=v2/.sol//First

E1=Grad[-v11,{R,θ,φ},"Spherical"]//

Simplify

E2=Grad[-v22,{R,θ,φ},"Spherical"]//

Simplify

E2=Grad[-v22,{x,y,z},"Cartesian"]

?VectorFieldPlots′;

GradientFieldPlot[-v3,{z,-2,2},{x,-2,2},PlotPoints→20]

Plot3D[v3,{z,-2,2},{x,-2,2},PlotPoints→30,LabelStyle→(FontSize→20),BoxRatios→{1,1,1}]

Mathematica繪制出的電場線圖和等電位圖如圖2所示.

(a)電場線圖

接下來對這題做一些極限情況下的討論：

(1)第一種極限情況：當ε=ε0時，介質球與環境的介電常數一樣，其對外場的響應消失，把其代入式(22)、(23)，此時φ1=φ2=-E0rcosθ空間的電場就是均勻電場.用Mathematica來分析此時的電場，繪制電場線圖的主要代碼如下所示，繪制的電場線圖如圖3所示.

Clear["Global′*"]

v11=-E0RCos[θ]

v22=-E0RCos[θ]

Grad[-v11,{R,θ,φ},"Spherical"]//Simplify

E2=Grad[-v22,{R,θ,φ},"Spherical"]//

Simplify

v22=-E0x;

E2=Grad[-v22,{x,y,z},"Cartesian"]

?VectorFieldPlots′;

GradientFieldPlot[-v3,{z,-2,2},{x,-2,2},PlotPoints→20]

圖3 極限情況1下介質球附近的電場線圖

結合φ1=φ2=-E0rcosθ和圖3進一步驗證了此種情況下的電場為勻強電場.

(24)

介質球的總電偶極矩為

(25)

式(22)φ1中的第二項正是這個電偶極矩所產生的電勢

(26)

由于此時的外電場相當于偶極子產生的電場，所以用Mathematica來繪制偶極子的電場來近似這種情況下的電場，Mathematica的代碼如下所示：

Clear["Global′*"];

V=q/(4*Pi*Subscript[ε,0])*(1/r1-1/r2)

cond={r1→R-l/2*Cos[θ],r2→R+l/2*

Cos[θ]}

V=V/.cond//SimplifySeries[V,{l,0,6}]//

Normal

Solve[V1==c,R]

Ee=Grad[-V1,{R,θ,ω},"Spherical"]

D[R[θ],θ]/R[θ]==Ee[[1]]/Ee[[2]]

DSolve[%,R[θ],θ]

g2=PolarPlot[Evaluate[Table

Identity];

Show[{g1,g2},DisplayFunction→$Display

Function,Epilog→{PointSize[0.03],Hue[2],

Point[{0,-2.5}],PointSize[0.03],Hue[.3],

Point[{0,2.5}]}]

V1[θ_,R_]:=cos[θ]/R2

RevolutionPlot3D[V1[θ,R],{R,0.01,3},{θ,0,2π},PlotRange→{-3,3},PlotPlints→30,

BoxRatios→{1,1,1}]

Mathematica繪制出的第二種極限情況下也就是電偶極子的電場線圖和等位面圖如圖4所示.

從圖4(a)中可以看出電偶極子的整個電場和電勢是一個數值分布對稱的區域，整個圖形是以電偶極子的軸線對稱的.再看圖4(b)，發現在正電荷附近電偶極子的電勢比較高；而在負電荷附近，電偶極子的電勢比較低，且在電偶極子的中垂線上電勢為零.

(a)電場線圖

(b)等位面圖

(3)當ε→-∞時，介質球內的場為E內→0，其效果相當于一個導體球.而此時，p=4πε0a3E0，問題就變成了求解導體球附近的電場，所以，我們可以用導體球的情況來近似這種情況下的電場分布，求導體球附近的電場的Mathematica代碼如下所示：

Clear["Global′*"]

v=Sum[(An*R^n+Bn/R^(n+1))*Pn,{n,0,∞}]

v=Sum[(An*R^n+Bn/R^(n+1))*Pn,{n,0,1}]/.{P0→1,P1→Cos[θ]}

(v/.{B0→0,B1→0})==-E0*R*Cos[θ]

cond1={A1→-E0,A0→0};

(v/.cond1/.R→a)==0

cond2={B1→a^3E0,B0→0};

v=v/.cond1/.cond2

GradientFieldPlot[-v1[z,x],{z,-2.,2.},{x,-2.2.},PlotPoints→20,Epilog→{Hue.

[.6],Disk[{0,0},1]}]

Mathematica繪制出的第三種極限情況下也就是導體球附近的電場線如圖5所示.

圖5 導體球附近的電場線圖

從圖中我們發現，導體球外的電場分布規律和介質球外的電場分布規律是一樣的，但是在這種情況下，導體球內部是沒有電場存在的.因為當導體置于電場中時，導體內的電荷會重新分布，形成電場來抵消外電場，最終就會導致導體內的總電場為零.

4 結論

本文基于Mathematica討論了分離變量法在靜電場的應用，通過用分離變量法求解介質球附近的電場以及3種極限情況下電場分布的實例分析展示了Mathematica具有的強大的數值計算、符號計算和圖形繪制的功能，研究探討了如何通過Mathematica軟件輔助靜電場的學習.學生通過學習Mathematica軟件，能夠很好地鍛煉動手能力和自主學習能力，他們在學習過程中可以將專業課的復雜數學問題交給Mathematica來進行，而自己只需掌握基本的數學原理，了解相關知識，配合Mathematica圖形顯示，就能達到更充分更深層次理解內容本質的目的，尤其是西藏地區來自農牧區的學生，數理基礎很差，Mathematica的引入解決了這些學生由于薄弱的數學基礎和抽象的專業概念所引起的在專業課學習上的困難，對這些學生具有重大意義.