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

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