. 第一類邊界問題的有限差分法探討 摘要：本次重點是對于第一類邊界問題的兩種不同方法的對比研討，通過計算機仿真有限差分法和計算分離變量法對同一問題的求解，對結(jié)果進行對比，能夠發(fā)現(xiàn)有限差分法更加快捷簡便，只要迭代次數(shù)足夠多就能使誤差趨于零。而分離變量法則是準確的計算出結(jié)果，只是運算相對復雜。 關(guān)鍵字:有限差分法，分離變量法，加速收斂因子，迭代次數(shù)，邊界條件。 引言：在給定的三類邊界條件①下求解標量位或矢量位的泊松方程或拉普拉斯方程的解一般的理論依據(jù)是唯一性定理和得加原理，由此而得出的解題方法有很多。主要分為兩大類：一是解析法（如分離變量法，鏡像法②等），二是數(shù)值法（如有限差分法，有限元法③等）。這兩種方法各有優(yōu)點和不足④，相比較而言在許多實際問題中由于邊界條件過于復雜而無法求得解析解。這就需要借助于數(shù)值法來求電磁場的數(shù)值解。有限差分法便是一種比較容易的數(shù)值解法。本次研討就以第一類邊界問題進行為例來分析研究有限差分法。 一、 有限差分法的定義： 微分方程和積分微分方程數(shù)值解的方法為有限差分法?；舅枷胧前堰B續(xù)的定解區(qū)域用有限個離散點構(gòu)成的網(wǎng)格來代替， 這些離散點稱作網(wǎng)格的節(jié)點；把連續(xù)定解區(qū)域上的連續(xù)變量的函數(shù)用在網(wǎng)格上定義的離散變量函數(shù)來近似；把原方程和定解條件中的微商用差商來近似， 積分用積分和來近似，于是原微分方程和定解條件就近似地代之以代數(shù)方程組，即有限差分方程組 ， 解此方程組就可以得到原問題在離散點上的近似解。然后再利用插值方法便可以從離散解得到定解問題在整個區(qū)域上的近似解。 精品 . 二、 有限差分法解題的基本步驟： （1）、區(qū)域離散化，即把所給偏微分方程的求解區(qū)域細分成由有限個格點組成的網(wǎng)格； （2）、近似替代，即采用有限差分公式替代每一個格點的導數(shù)； （3）、逼近求解。換而言之，這一過程可以看作是用一個插值多項式及其微分來代替偏微分方程的解的過程。 三、 有限差分法公式的推導： 把求解的區(qū)域劃分成網(wǎng)格，把求解區(qū)域內(nèi)連續(xù)的場分布用網(wǎng)絡節(jié)點上的離散的數(shù)值解來代替。網(wǎng)格劃分的充分細，才能夠達到足夠的精度。應用有限差分法計算靜態(tài)場邊值問題時，需要把微分方程用差分方程替代。 用圖形法解釋如下： 精品 . 有限差分網(wǎng)格劃分 在由邊界L界定的二維區(qū)域D內(nèi)，電位函數(shù)φ滿足拉普拉斯方程且給定第一邊界條件，則： 如圖將區(qū)域D劃分為正方形網(wǎng)格，網(wǎng)格線的交點稱為節(jié)點，兩相鄰平行網(wǎng)格線間的距離稱為步距 h。然后，拉普拉斯方程離散化，對于任一點0，有一階偏導數(shù)： 而后，對于二階偏導數(shù)： 精品 . 對于Y軸同理： 因此拉普拉斯方程的差分格式為： 緊鄰邊界節(jié)點的拉普拉斯方程的差分格式為： 其中p、q為小于1的正數(shù)；1、2為邊界上的節(jié)點，其值為對應邊界點處的值，是已知的。具體如圖： 緊鄰邊界節(jié)點的網(wǎng)格劃分 精品 . 應用數(shù)值計算解釋（泰勒公式展開法）： 1點電位的泰勒公式展開為 3點電位的泰勒公式展開為 ，當h很小時，忽略4階以上的高次項，得 同樣可得 將上面兩式相加得 在上式中代入，得 精品 . 對于，即F=0的區(qū)域，得到二維拉普拉斯方程的有限差分形式 通過以上兩種方法的推導，可得任意點的電位等于圍繞它的四個點的電位的平均值。當用網(wǎng)格將區(qū)域劃分后，對每一個網(wǎng)絡點寫出類似的式子，就得到方程數(shù)與未知電位的網(wǎng)絡點數(shù)相等的線性方程組。已知的邊界條件在離散化后成為邊界上節(jié)點的已知電位值。 四、 差分方程組的解法 方法一：高斯——賽德爾迭代法（簡單迭代法） 其步驟是先對每一網(wǎng)格點設一初值。然后按一個固定順序（一般點的順序按“自然順序”，即：從左到右，從下到上）如圖： 之后，利用二維拉普拉斯方程的有限差分形式用圍繞它的四個點的電位的平均值作為它的新值，當所有的點計算完后，用它們的新值代替舊值，即完成了一次迭代。然后再進行下一次迭代，如此循環(huán)。如下式： 精品 . （迭代公式1） 其式中的上角標（k）表示k次近似值，下腳標i，j表示節(jié)點所在位置，即第i行第j列的交點。其中要特別注意：在迭代過程中遇到邊界點式，需將邊界條件 帶入。 循環(huán)迭代時當所有內(nèi)節(jié)點滿足以下條件時停止迭代： 其中，W是預定的最大允許誤差。 方法二：逐次超松弛法 簡單迭代法在解決問題時收斂速度比較慢，實用價值不大。實際中常采用逐次超松弛法（又稱高斯——賽德爾迭代法變形），相比之下它有兩點重大的改進 ，第一是計算每一網(wǎng)格點時，把剛才計算得到的臨近點的新值代入，即在計算(i,j)點的電位時，把它左邊的點(i-1,j)和下面的點(i,j-1)的電位用剛才算過的新值代入，即: （迭代公式2） 第二，是引入“加速收斂因子”。上式中的α即為“加速收斂因子”，且1《 α <2。特別關(guān)注的是逐次超松弛法收斂的快慢與α有明顯關(guān)系。并且最佳α的取值隨著條件的不同而不同，如何選擇最佳α，是個復雜問題。 精品 . 在計算時可以嘗試求取最佳α值，以使計算快速。 五、 應用計算機仿真有限差分法解決具體問題 本次討論我選擇第四題為具體實例進行研究。 題目： 如圖所示，有一長方形的導體槽，a = 20，b = 15，設槽的長度為無限長，槽上有一塊與槽絕緣的蓋板，電位為100V，其他板電位為零，求槽內(nèi)的電位分布。 通過MATLAB進行仿真，運用有限差分法，源代碼如下： u=zeros(15,20); i=2:14; u(i,20)=100; for j=1:20 for i=2:14 u(i,j)=100/19*(j-1); end end a=input(please input a(1<a<2);a=); m=input(please input m(1<m);m=); 精品 . for n=1:m; for i=2:14; for j=2:19; b=u(i,j); c=u(i,j+1); d=u(i+1,j); e=u(i-1,j); f=u(i,j-1); g=0.25.*(c+d+e+f); u(i,j)=b+a.*(g-b); end end end mesh(u); 通過對相關(guān)資料的查詢及學習，我認為運用MATLAB編程就是要將有限差分法的基本公式實現(xiàn)。以此為出發(fā)點可以讓解析和編寫過程相對簡單化。 首先，先要對電場的邊界條件進行設定： u=zeros(15,20); i=2:14; u(i,20)=100; 使其滿足題目的要求右側(cè)邊界電位為100，其余三邊為0，由此也可以判斷出此題目為第一類邊界條件的問題。 精品 . 同時，此處，也要特別注意，因為在MATLAB程序中，編寫一個矩陣的默認順序是從左到右，從上到下。另外，命名矩陣時是先編寫行，后編寫列，因此，對于u(i,j)中，i代表的是第i行即對應縱坐標，j代表的是第j列即對應橫坐標，與坐標表示有一定的差別，需要區(qū)分和注意。 之后是相對關(guān)鍵的編寫部分，就是對電場內(nèi)部的電位設定，此處，我通過分離變量法將u分解為u（x）和u(y)，分別求解。為讓計算簡便，我設定其內(nèi)電場為線性變化的，即u（x）=kx+b，u（y）=0作為初值進行計算（應用分離變量法得到的結(jié)果是u(x)=Bsinh m∏/15 u(y)=Asinm∏/15。可以通過此處在得出結(jié)果后進行對比和誤差分析。)又因為要在MATLAB中進行運算，所以首先要轉(zhuǎn)化為i和j的形式，因此可得：u(j)=kj+b和u(i)=0。從而得到方程組： u(1)=k+b=0; u(20)=20k+b=100。 解得：u(j)=100/19*(j-1),即：u(i,j)=100/19*(j-1)。由此得到了內(nèi)電場除去上下邊界（衡為零）的電位方程，在MATLAB中我運用了for循環(huán)語句來實現(xiàn)其內(nèi)電場的賦值，源代碼如下： for j=1:20 for i=2:14 u(i,j)=100/19*(j-1); end 精品 . end 將電場各點的電位賦初值完成后，就需要帶入公式進行計算了。在之前先要對“加速收斂因子” α的最佳值和預定的最大允許誤差進行確定。關(guān)于加速因子的最佳值確定問題，我通過學習和查閱相關(guān)資料得到了兩種方法，如下： 一是逐步搜索法； 將α的取值區(qū)間（1,2）進行M等分，α分別取1+1/M,1+2/M……1+（M-1）/M,通過上面提到的迭代公式2依次對同一個精度要求求出迭代次數(shù)k的值，并從中選出“加速收斂因子” α的最佳值，具體步驟如下： 步驟1 給定等分數(shù)M和精度要求ε的值，令α的初始值為1； 步驟2 令p=1,2,3,……M-1，重復步驟3-5； 步驟3 αp=1+p/M; 步驟4 按照如下公式迭代 找出符合精度要求ε的迭代次數(shù)kp; 步驟5 比較找出kp值最小的αp為最佳的加速收斂因子α的值。 二是黃金分割法： 依據(jù)黃金分割法的思想，通過計算機自動選取最佳加速收斂因子α的近似值，具體步驟如下： 步驟1 對（1,2）區(qū)間第一次0.618的分割，區(qū)間分界a 精品 . 1=1,b1=2,在（a1，b1）區(qū)間分割出黃金點point1= a1+0.618(b1- a1),進行迭代公式2的迭代，求出迭代次數(shù)K值，如果迭代次數(shù)沒有超出規(guī)定的發(fā)散常數(shù)，迭代結(jié)束，否則轉(zhuǎn)為步驟2。 步驟2 在（1,1.618）和（1.618,2）之間進行第二次的黃金分割，找出分割點point2=a2+0.618(b2- a2),其中a2和b2是新分割區(qū)間的左右邊界。找出迭代次數(shù)最少的α。 以上兩種方法是比較具有實踐和研究價值進行的，在這里介紹和分享給大家探討。在我分析過后認為第二種較好，因為其可以通過計算機自動選取到最佳的加速收斂因子，但是實施起來比較復雜，要通過MATLAB或其它語言進行編碼實驗得到，我初步試驗過，沒能成功，就運用了第一種方法進行了計算，但是，我發(fā)現(xiàn)其運算過于復雜，筆算正確的成功率比較低，而運用電腦編程要重復多次運算，雖然，計算機運算快速且方便，但是，這樣的源代碼同樣過于復雜，可以進一步進行深入的研究。 本次的題目解答中，由于是應用計算機計算，所以計算速度和計算效率在不同的α取值中沒有明顯區(qū)別。因此，我就只是設定了一個1~2之間的一個數(shù)字：1.8為“加速收斂因子” α的值了，進行了近似計算，實驗誤差在可接受的范圍。 然后，是對于迭代次數(shù)的確定。當?shù)鷿M足停止條件： 時，對應的迭代次數(shù)即為最終的迭代次數(shù)，這里的W為設定好的允許誤差，為從簡方式，我選擇運用輸入迭代次數(shù)，并通過多次輸入對比輸出結(jié)果選取最佳迭代次數(shù)的方法來選取，存在了一定的誤差。但是，相比于我運用上面的公式求解計算迭代次數(shù)，進而得出的圖形的誤差比要小許多，我一直嘗試進行改進，但一直沒能很好的解決誤差過大的問題：最主要的問題是邊界值一直不滿足我的設定值。希望通過進一步的研討和大家的交流能過得到最優(yōu)方式。 精品 . 通過以上的分析，我得出了進一步的編程代碼： a=input(please input a(1<a<2);a=); m=input(please input m(1<m);m=); 其中a和m分別為最佳“加速收斂因子”和迭代次數(shù)，均為我們輸入數(shù)值計算的方式。先要預算出一些取值范圍，然后通過分別輸入，對比輸出結(jié)果進而來選取最優(yōu)解。 之后的源代碼的編寫思想就是通過迭代計算公式來實現(xiàn)對電場內(nèi)部各點電位值的計算，因此，在輸入“加速收斂因子”和“迭代次數(shù)”的數(shù)值后，我通過for的循環(huán)語句實現(xiàn)了對電場內(nèi)部的電位迭代求解，并用圖形形式表示出來。其源代碼如下： for n=1:m; for i=2:14; for j=2:19; b=u(i,j); c=u(i,j+1); d=u(i+1,j); 精品 . e=u(i-1,j); f=u(i,j-1); g=0.25.*(c+d+e+f); u(i,j)=b+a.*(g-b); end end end mesh(u); 上述代碼中，b,c,d,e,f分別場內(nèi)一點及其左右上下各點的電位值，之后的g和u（i，j）則實現(xiàn)了逐次超松弛法的迭代，最后用mesh實現(xiàn)圖形的輸出。 六、 結(jié)果和誤差的對比分析 由分離變量法得到的初步結(jié)果是u(x)=Bsinh m∏/15和u(y)=Asinm∏/15，從而應有u=Csinh(m∏/15)sinm(∏/15)的結(jié)果，從其表達式的形式，我們可以簡單推斷出電場的電位在x軸上的分布是成虛指數(shù)函數(shù)的形式呈現(xiàn)，而在y軸上是成正弦函數(shù)呈現(xiàn)的，這和通過有限差分法的計算機仿真的圖形基本符合，如下圖： 精品 . 因此，可以基本確定通過計算機仿真得到的結(jié)果的正確性。 另外，可能存在較大誤差的方面是跌代次數(shù)，因為迭代次數(shù)是我們?nèi)斯ぽ斎氲?，所以判斷其正確與否要通過輸出的結(jié)果來判定、來選擇，下面幾幅圖是當“加速收斂因子” α=1.8時，不同的迭代次數(shù)m的對應圖形： m=1時： m=10時： 精品 . m=20 m=30 m=32 精品 . m=35 m=40 通過上面迭代次數(shù)m取不同的值對應的不同圖形分析可以看出，在m=1,10,20時期計算結(jié)果都存在較大誤差，圖形與分離變量法得到的結(jié)果相差很大，有明顯的不規(guī)則部分，不成現(xiàn)光滑的曲線。而當 精品 . m=30后圖形與分離變量法計算到的圖形近似，而且圖形規(guī)則，曲線光滑，而且隨著m值的增大圖形的變化非常微小了，因此，我選擇m=32為迭代次數(shù)，所得到的結(jié)果誤差比較小，結(jié)果相對準確。 七、 研究結(jié)論 學習了運用有限差分法計算第一類邊界條件的電位值，通過計算機仿真進行了實際例題的求解，通過與分離變量法求解的對比，證明可行性和有效性，適用于解決第一類邊界問題，而且更加方便和快捷。同時，“加速收斂因子”最佳解問題比較復雜，不過通過計算機仿真的過程中，不同數(shù)值“加速收斂因子”的加速效果不易體現(xiàn)出來。跌代次數(shù)對計算結(jié)果又較大影響，迭代次數(shù)越多，結(jié)果越準確，誤差越小，但相應的計算也會更加復雜，因此只要迭代次數(shù)得到得結(jié)果滿足一定的允許誤差范圍即可。 附錄： ① 三類邊界條件：第一類：已知電位在邊界上的數(shù)值。 第二類：已知電位的導數(shù)在邊界上的數(shù)值。 第三類：已知電位在一部分邊界上的數(shù)值和在其余邊界上的導數(shù)的數(shù)值。 ② 解析法包括：分離變量法，鏡像法，復變函數(shù)法，保角變換法，格林函數(shù)法。 ③ 數(shù)值法包括：數(shù)值積分法，有限差分法，有限元法，矩量法。 ④ 解析法的優(yōu)點：由有限個項構(gòu)成閉合解或無窮級數(shù)，通常具有鮮明的物理意義。 精品 . 缺點：解題范圍窄小，對邊界形狀十分挑剔。 數(shù)值法的優(yōu)點：解題范圍寬廣，對邊界的形狀沒有限制。 缺點：數(shù)據(jù)離散，物理意義深藏其中。 心得體會： 有限差分法的公式套用很簡單，記住結(jié)果很簡單，但是要想深刻理解其內(nèi)容和靈活應用則比較困難，需要一定時間的深刻學習和做題應用，同時，這個方法確實是一個非常方便和有效地解決三類邊界條件問題的工具。 MATLAB在電磁場方面的應用相當廣泛，其內(nèi)擁有的很多函數(shù)都可以模擬靜電場中的實例，就比如這次的研討就應用了矩陣，通過有限差分法來模擬靜電場的電位分布，相信將來還會有更多的有意思，有研究價值計算和模擬電磁場的相關(guān)知識。是學好電磁場必須要牢固掌握的工具。 如有侵權(quán)請聯(lián)系告知刪除，感謝你們的配合！ 精品