基于Mathematica和Matlab的典型數(shù)學(xué)物理方程解分布圖像制作.doc

第 29 卷 第 7 期2012 年 7 月吉林化工學(xué)院學(xué)報Vol 29 No 7JOURNAL OF JILIN INSTITUTE OF CHEMICAL TECHNOLOGYJul2012文章編號:1007-2853( 2012) 07-0041-04基于和的典型數(shù)學(xué)物理MathematicaMatlab方程解分布圖像制作陳殿偉1 ，楊海英2( 1 吉林化工學(xué)院 理學(xué)院，吉林 吉林 132022; 2 吉林市第五中學(xué) 物理組，吉林 吉林 132011)摘要: 通過計算機(jī)輔助軟件 Mathematica 和 Matlab，把典型數(shù)學(xué)物理方程解的空間分布制作成三維圖像關(guān) 鍵 詞: Mathematica; Matlab; 數(shù)學(xué)物理方程中圖分類號: O 411 1文獻(xiàn)標(biāo)志碼: A典型數(shù)學(xué)物理方程包括: 波動方程; 熱傳導(dǎo)方程; 泊松方程，u = 0 為拉普拉斯方程，二者都是 穩(wěn)定分布方程 這三個方程構(gòu)成了“數(shù)理方程”的 主要內(nèi)容 它們都是二階線性方程，刻畫了很多物 理現(xiàn)象的規(guī)律1-3波動方程分布圖像1波動方程或稱波方程是一種重要的偏微分方程，主要描述自然界中的各種的波動現(xiàn)象 其表達(dá) 式為:圖 1 位移在時間和空間中的分布圖像2( 1)utt = a u + f 熱傳導(dǎo)方程2如弦的 x = 0 端固定，x = l 端受迫作諧振動 F= Asint，弦的初始速度為零，求弦的振動 這個 定解問題是:當(dāng)一個物體內(nèi)部各點(diǎn)的溫度不一樣時，則熱量就會從溫度高的地方向溫度低的地方流動，這2utt a uxx = 0 ，u | x = 0 = 0，u | x = l = Asint ，u | t = 0 = 0，ut | t = 0 = A x = l 端為非齊次邊界條件種現(xiàn)象就是熱傳導(dǎo)問題 由于熱傳導(dǎo)過程總是表( 2)現(xiàn)為溫度隨時間和位置的變化，所以，解決熱傳導(dǎo)問題歸結(jié)為求物體內(nèi)部溫度的分布問題 如果研究物質(zhì)的擴(kuò)散問題，物質(zhì)的擴(kuò)散是由濃度高的地這是一個不受外力作用的振動，但它 在 弦x = l的一端，有一個諧振源 F = Asint，而在初始方擴(kuò)散到低的地方，且服從傅立葉熱傳導(dǎo)方程相類似的能斯特擴(kuò)散定理，即物質(zhì)的擴(kuò)散與濃度的時刻，速 度 為 A 用行4，得到如下圖形:Mathematica直 接 編 程 運(yùn)變化成正比 與熱傳導(dǎo)方程類似，可以得到濃度 u滿足熱傳導(dǎo)方程，所以熱傳導(dǎo)方程也稱擴(kuò)散方程其表達(dá)式為圖 1 中，a = 2，A = 5，l = 1， = 1，x 0，1 ，t0，10 從圖中看出，在弦 x = 0 的一端，由于沒 有諧振源，因此，位移為 0; 而在 x = l 的一端，諧振 源 F = Asint，位移 u( x，t) 隨著 F 的作用而變化= a2 u u( 3)t對于細(xì)桿導(dǎo)熱問題 初始時刻桿的一端溫度為零度，另端溫度為 u ，桿上溫度梯度均勻，零度0收稿日期:2012-04-14作者簡介:陳殿偉( 1968-) ，男，吉林省吉林市人，吉林化工學(xué)院副教授，主要從事凝聚態(tài)物理方面的研究的一端保持溫度不變，另一端跟外界絕熱桿上溫度 u ( x，t) 滿足下列泛定方程和定解 條件化，因?yàn)闂U上的每個領(lǐng)域溫度已經(jīng)達(dá)到平衡穩(wěn)定分布方程分布圖像32ut a uxx = 0 ，u | x = 0 = 0，ux | x = l = 0 ，在熱傳導(dǎo)問題中，如果溫度分布穩(wěn)定，熱源強(qiáng)度 f( x，y，z) 不隨時間變化，熱傳導(dǎo)持續(xù)下去，最終 將達(dá)到穩(wěn)定狀態(tài)，空間中各點(diǎn)的溫度不再隨時間( 4)2u | t = 0 = u0 x / l ( a= k / c)泛定方程和邊界條件都是齊次的，用分離變數(shù)法 求 解 在 這 里 可 以 用 Mathematica 直 接 編 程4，就可以畫出桿上溫度 u ( x，t) 的分布圖像， 如圖 2 3 所示3變化 ，即 ut = 0，得到方程u + f( x，y，z)= 0 ，( 5)此方程稱為泊松方程，如果沒有熱源，即 f = 0則得到拉普拉斯方程 在擴(kuò)散問題中，濃度處于穩(wěn) 定的狀態(tài)考慮振動的平衡現(xiàn)象，同樣得到穩(wěn)定分 布方程3 1 拉普拉斯方程如帶電的云跟大地之間的靜電場近似是勻強(qiáng) 靜電場，其電場強(qiáng)度 E0 是豎直的 水平架設(shè)的輸 電線處在這個靜電場之中，如圖 4 所示 輸電線是 導(dǎo)體圓柱，柱面由于靜電感應(yīng)出現(xiàn)感應(yīng)電荷，圓柱 鄰近的靜電場也就不再是勻強(qiáng)的了 不過，離圓柱 “無限遠(yuǎn)”處的靜電場仍保持為勻強(qiáng)的 現(xiàn)在研究 導(dǎo)體圓柱怎樣改變了勻強(qiáng)靜電場圖 2 溫度在時間和空間中的分布圖 3 溫度在時間和空間中的分布圖 2 3 的視角不同，x 0，1 ，t 0，5 ，從兩圖形中可以看出桿上 x = 0 的一端，由于受其 邊界條件的限制，u( x，t) 不隨時間的變化而變化， 它仍然是零度; 在 x = l 的一端，t = 0 時 u( x，t) 在最高點(diǎn)，隨著時間的增大，u( x，t) 迅速減小，變化 的幅度很大，然后趨于平緩，最后幾乎沒有什么變 化，但 t = 0 時，u ( x，t ) 是從桿上 x = l 的一端向 x = 0的一端傳導(dǎo)的，只是變化幅度不是很大 從整 個圖形來看，u( x，t) 除了在邊界 x = 0 處，它都是 隨著時間的增大而迅速減小，然后趨于平緩，到后 來根本沒有什么變化 這很明顯是一個熱傳導(dǎo)方 程的圖形，溫度從高點(diǎn)向低點(diǎn)飛快的傳導(dǎo)，隨著時 間的增長，桿上的溫度越接近均勻，所以我們后來圖 4 電場強(qiáng)度在空間的分布首先需要把這個物理問題表明為定解問題取圓柱的軸為 z 袖 如果圓柱“無限長”，那么，這 個靜電場的電場強(qiáng)度、電勢顯然跟 z 無關(guān)，我們只 需在 oxy 平面上加以研究就夠了如圖 4 所示，正是 oxy 平面上的靜電場，圓柱 面在 oxy 平面的剖口是圓 x2 + y2 = a2 ，其中 a 是 圓柱的半徑柱外的空間中沒有電荷，所以電勢 u 滿足二 維的拉普拉斯方程:uxx + uyy = 0( 柱外)( 6)導(dǎo)體中的電荷既然不再移動，這說明導(dǎo)體中第 7 期陳殿偉，等: 基于 Mathematica 和 Matlab 的典型數(shù)學(xué)物理方程解分布圖像制作43全可以把導(dǎo)體的電勢當(dāng)作零，從而邊界條件為面均勻地一層層切割 u 所得到的，所以，由等高線u | 2 + y2 = a2 = 0 ( 7)組成的圖形就像是 u 在 oxy 平面上的倒影x按照分離變量法 u( x，y) = X( x) Y( y) 代入拉普拉斯方程固然不難把它分解為兩個常微分方 程，但代入上述邊界條件卻只能得到:X( x) Y( 槡a2 x2 ) ，( 8)不能分解為 X( x) 或 Y( y) 的邊界條件 事實(shí)上，既然邊界是圓，直角坐標(biāo)系顯然是不適當(dāng)?shù)?，必須采用平面極坐標(biāo)系拉普拉斯方程在極坐標(biāo)系中的表達(dá)式為2 u / 2 + ( 1 / ) u / + ( 1 / 2 ) 2 u / 2 = 0，圖 5 極坐標(biāo)中的電勢分布圖像( a)，( 9)由( 9) 式推得其電場強(qiáng)度的表達(dá)式為式中: 是極徑， 是極角 “導(dǎo)體電勢為零”就表明為齊次的邊界條件:E = ( u cos usin)+ ( u sin + ( 10)u | = a = 0 u cos，在“無限遠(yuǎn)”處的靜電場仍然保持為勻強(qiáng)的E0 由于選取了 x 鈾平行于 E0 ，所以在無限遠(yuǎn)處，Ey = 0，Ex = E0 ，即 u / x = E0 ，亦即 = ( D2cosx + Ea 22 E2 ( sin x cos x) ) +000yysinx E0 cos 因而還有一個非齊次u = E0 x =的邊界條件:( D) ( 13)0y由( 13) 可求出 E，利用 Matlab 編程，并繪制( 11)u | = E0 cos 求得柱外的靜電勢為:圖形，如圖 6 所示2u( ，) = D0 ln ( / a) E0 cos E0 ( a / )cos ，( 12)其中: E0 cos 項為原勻強(qiáng)靜電場中的電2勢分布 當(dāng) 較大時， E0 ( a / ) cos 可以忽略，所以它是對圓柱附近勻強(qiáng)電場的修正項，原因是 受柱面感應(yīng)電荷的影響 此外，還有 D0 ln ( / a ) 項，它的系數(shù) D0 是任意常數(shù)，這說明包含著某種 不確定的因素，這個不確定因素在于問題提出時 根本沒有說明導(dǎo)體柱原來所帶的電 量，可 見 D0 ln( / a) 正是圓柱原來所帶電量的影響( 由靜 電學(xué)可知，D0 ln( / a) 正是均勻帶電圓柱體周圍的 靜電場中的電勢) 根據(jù)其最終解( 12) ，用 Matlab 直接編程并運(yùn) 行5-6，如圖 5 所示圖 5 中的 D0 = 3，a = 1，E0 = 5，0，2， =1，10，這個圖就像是一個圓柱上面頂著一 個圓面，柱和面接觸的地方不動時，圓面整個彎曲 了一些，就形成上面的這個圖 從圖中看到其邊界 是個圓，正好符合了題中給出的邊界條件是個圓 的要求 從上往下看，電勢分布的方向由高向低 的，但是，當(dāng)電勢降低到圓中心處時，圓的中心是空 的，因?yàn)閳D中圓的中心處正好是導(dǎo)體圓柱所在的地 方，導(dǎo)體內(nèi)部是個等勢體，故變化是零 圖下面的線 條是電勢的等高線，這些等高線是由平行于 oxy 平圖 6 極坐標(biāo)中的電場分布圖 6 是關(guān)于電場強(qiáng)度方向的圖形，其中 D0 = 3a = 1，E0 = 5，0，2，1，10 ，我們看到的這個圖形是在極坐標(biāo)下表示的圖形，在這個圖中，左邊電力線從無窮遠(yuǎn)處向近處平行的行進(jìn)著可是到達(dá)中心時，沒有了電力線 這是因?yàn)閷?dǎo)體圓柱體正在這個中心上，由于導(dǎo)體內(nèi)部為等勢體，所以電場強(qiáng)度為零，也就沒有了電力線，圖上的電力線從遠(yuǎn)到近方向行進(jìn)，結(jié)束于導(dǎo)體的表面，后又從導(dǎo)體的表面出發(fā)，向遠(yuǎn)處行進(jìn)泊松方程在圓域 0 上求解泊松方程的邊值問題3 2u = a + b( x2 y2 ) ，u |根據(jù)其最終解:= c ( 14) = 0u = + = c + a ( 2 2 ) + b 2 ( 2 2 )也有兩個相等的最小值，而電勢是沿著電荷密度分布的，電荷密度越大的地方，其電勢也越大，所 以，電勢分布圖中，有兩個相等的最大值，和兩個 相等的最小值 沿著 y 軸，有兩個最大值，沿著 x 軸有兩個最小值 由于邊界條件的限制，其電勢分 布邊界為一個圓，而在這個圓的平面上，電勢為 零 下面那些是等高線，它們是平行 oxy 平面切割 得到的00412( 15)cos2 利用 Matlab 直接編程5-6，就可以把函數(shù) u繪成一個三維圖形，如圖 7 8 所示結(jié)論4主要針對數(shù)學(xué)物理方程中典型的數(shù)學(xué)物理方程，使用 Mathematica 和 Matlab 軟件進(jìn)行編程繪 制出位移的分布圖像; 對齊次和非齊次的熱傳導(dǎo) 方程，利用了中的 PDE 工具箱，通過相應(yīng)的數(shù)值 計算步驟，畫出了溫度的分布圖像，對于穩(wěn)定分布 方程，制作出了物理量( 標(biāo)量) 和相應(yīng)向量場的分 布圖像圖 7 電勢分布圖像參考文獻(xiàn):1梁昆淼 數(shù)學(xué)物理方法M 北京: 高等教育出版社，2010姚端正，梁家寶 數(shù)學(xué)物理方法M 北京: 人民教 育出版社，1997黃大奎，舒慕曾 數(shù)學(xué)物理方法M 北京: 高等教 育出版社，2001吳劍，胡波主 掌握 和 精 通 Mathematica4 0M北京: 北京郵電出版社，2002陸 君 安，等 偏 微 分 方 程 的 MATLAB 解 法M武漢: 武漢大學(xué)出版社，2001劉宏友，李莉，彭鋒 MATLAB 6 基礎(chǔ)及應(yīng)用M重慶: 重慶大學(xué)出版社，20012圖 8 電荷密度分布圖像3圖 7中，c= 3，a = 1，b= 4，0= 5， 80，2，1，10 ，這兩個圖像都是由極坐標(biāo)轉(zhuǎn)換成直角坐標(biāo)，第一個是電勢分布問題，而第 二個圖像是為了更好地說明第一個圖形而添加的 電荷密度分布問題，電荷密度 w = a + b ( x2 y2 ) 轉(zhuǎn)換成極坐標(biāo)為 w = a + b( 2 cos2) ，從電荷密度 的公式及其圖像，可以看出 w 有兩個相等最大值，456Fabrication of Distribution Image for Solution of Typical MathematicalPhysics Equation Based on Mathematica and MatlabCHEN Dian-wei1 ，YANG Hai-ying2( 1 College of Sciences，Jilin Institute of Chemical Technology，Jilin City 132022，China; 2 Physics Group，NO 1 Mid