計(jì)算電磁學(xué)-第4章-有限差分法

第四章有限差分法差分與差商歐拉近似、梯形法則和龍格-庫塔法Poisson方程差分格式媒質(zhì)邊界的差分格式邊界條件的差分格式差分方程求解差分實(shí)例分析Helmholtz方程的差分法計(jì)算電磁學(xué)的應(yīng)用計(jì)算電磁學(xué)的應(yīng)用軍事方面的應(yīng)用計(jì)算電磁學(xué)的應(yīng)用軍事方面的應(yīng)用計(jì)算電磁學(xué)的應(yīng)用高速電子電路方面的應(yīng)用

計(jì)算電磁學(xué)的應(yīng)用超高速光集成電路的應(yīng)用計(jì)算電磁學(xué)的應(yīng)用超高速光集成電路的應(yīng)用cem@計(jì)算電磁學(xué)的應(yīng)用左手材料的仿真由蘇聯(lián)理論物理學(xué)家Veselago于1968年最早提出[1]，并預(yù)言了這種材料可能具有奇特的性質(zhì)。雙負(fù)材料的研究工作成為科學(xué)界關(guān)注的焦點(diǎn)，被美國(guó)Science雜志評(píng)為2003年度十大科技進(jìn)步計(jì)算電磁學(xué)的應(yīng)用左手材料的仿真計(jì)算電磁學(xué)的應(yīng)用隱身衣的仿真隱身斗篷的運(yùn)作秘訣就在于它能令微波的路徑變彎。這種斗篷其實(shí)以數(shù)千塊細(xì)小的“特異材料”片制成。這種人造纖維玻璃般的物料能控制光線。計(jì)算電磁學(xué)的應(yīng)用隱身衣的仿真cem@電磁仿真軟件介紹粒子仿真軟件：MAGICMAFIAFIT:有限積分法光學(xué)仿真軟件：ZEMAXGRASPWirelessInsitecem@cem@cem@cem@cem@cem@CST粒子仿真PierceGunＭＡＧＩＣcem@ＭＡＧＩＣ有限差分法介紹

有限差分方法是一種微分方法，自上世紀(jì)五十年代以來得到了廣泛應(yīng)用。概念清晰，方法簡(jiǎn)單，直觀。其與變分法相結(jié)合所形成的有限元法更有效。但有限差分還是以其固有特點(diǎn)在數(shù)值計(jì)算中有其重要地位，是應(yīng)用最多的一種數(shù)值方法。為求解由偏微分方程定解問題所構(gòu)造的數(shù)學(xué)模型，有限差分法是將定解區(qū)域（場(chǎng)區(qū)）離散化為網(wǎng)格離散節(jié)點(diǎn)的集合。以各離散點(diǎn)上函數(shù)的差商來近似該點(diǎn)的偏導(dǎo)數(shù)，使待求的偏微分方程定解問題轉(zhuǎn)化為一組相應(yīng)的差分方程。根據(jù)差分方程組解出各離散點(diǎn)處的待求函數(shù)值——離散解。

1、差分與差商

用差分代替微分，是有限差分法的基本出發(fā)點(diǎn)。這一點(diǎn)由微分原理保證的，當(dāng)自變量的差分趨于零時(shí)，差分變成微分

差分與差商

前向差分后向差分中心差分

差分與差商

通過泰勒公式分析上面差分精度，在某點(diǎn)上的一階導(dǎo)數(shù)的逼近度可由泰勒公式展開

兩式相減差分與差商

前向、后向差分截?cái)嘤?，具有h的一階精度，而中心差分法截?cái)嘤?，具有h的二階精度，中心差分的精度比較高。二階導(dǎo)數(shù)的差分格式前向差分后向差分都采用中心差分？差分與差商

對(duì)偏導(dǎo)數(shù)，可仿照上述方法，將表示為：2、歐拉近似、梯形法則和龍格-庫塔法

除了用最簡(jiǎn)潔差分去代替微分外，還可以從微分方程的解為積分形式出發(fā)，利用數(shù)值積分計(jì)算的手段處理差分格式，通過直觀圖形法去分析。設(shè)有初值問題歐拉近似

歐拉近似法在函數(shù)圖上用階梯的折線代替曲線

f(n)f(n+1)歐拉近似法的精度利用泰勒展開歐拉近似法的誤差量級(jí)為,為一階精度。歐拉近似法適用于區(qū)間不大、精度要求不高的場(chǎng)合。梯形近似

改善數(shù)值分析的精度關(guān)鍵是如何更好地近似曲邊梯形的面積。用函數(shù)曲線的端點(diǎn)連線代替曲線，而不是用矩形的平行線上式為梯形近似法，又稱改進(jìn)的歐拉折線法，改進(jìn)后的誤差為，具有二階精度，比歐拉折線法小一個(gè)數(shù)量級(jí)

（證明見P112）梯形近似

可以看出方程兩邊都包含n＋1步的待求未知量，是隱式格式，如何求解？先用歐拉折線法計(jì)算出0到n＋1步的未知值的預(yù)估計(jì)值，然后，再用梯形近似法計(jì)算。就是“預(yù)估計(jì)-校正法”。該方法截?cái)嗾`差仍然為O（h3）。

例題改進(jìn)歐拉法的近似計(jì)算由于求x=0.01的近似值解：法1：法2：龍格-庫塔法一般梯形近似法不能一步計(jì)算得到所要求的精度，需要采用多次迭代才能夠得到積分精度要求當(dāng)?shù)揭欢ǖ拇螖?shù)時(shí)，即滿足精度時(shí)，就認(rèn)為獲得所要求的精度積分。

龍格-庫塔法采用預(yù)估計(jì)-校正法的梯形近似法，只計(jì)算一次就能得到的截?cái)嗾`差為O（h3）結(jié)果，但當(dāng)期望一步到達(dá)r階精度，就發(fā)展了龍格-庫塔法

設(shè)K作為yn+1到y(tǒng)n的增量，得到n+1步長(zhǎng)處的結(jié)果K是用ω常數(shù)對(duì)K1,K2,…進(jìn)行線性組合構(gòu)造函數(shù)龍格-庫塔法K1,K2,…用α,β常數(shù)及K構(gòu)造的自變量表示的函數(shù)f的一系列值選取與函數(shù)f無關(guān)又與步長(zhǎng)數(shù)n無關(guān)的常數(shù)這些與計(jì)算函數(shù)無關(guān)的常數(shù)如何去確定？這也是龍格-庫塔法的關(guān)鍵，思路是：即要計(jì)算函數(shù)的增量，就要想到泰勒展開方法，上式中每個(gè)K值就是在展開處附近的近似值。

龍格-庫塔法右端(xn,yn)處做泰勒展開,重新按h冪整理，合并同類項(xiàng)，得

同樣對(duì)微分方程的解y(x)在點(diǎn)(xn,yn)進(jìn)行泰勒展開比較上面兩式，只要它們前面項(xiàng)的系數(shù)盡可能多的相等，就保證了截?cái)嗑取?/p>

龍格-庫塔法選取α、β、ω系數(shù)，使兩式項(xiàng)的系數(shù)相等

如果該關(guān)系式能夠一直維持到第m階仍能成立，但m＋1階不再成立，就稱為m階龍格-庫塔法二階龍格-庫塔法：二階龍格-庫塔法的形式解展開K2：則：則：將y(xn+1)進(jìn)行泰勒基數(shù)展開：比較以上兩式得：帶入上式，整理得龍格-庫塔法二階龍格-庫塔法的系數(shù)方程組這是一個(gè)不定方程組，也就是方程的變量的數(shù)目多于方程個(gè)數(shù)，因而方程的解有多個(gè)。上式如果取α2=β21=1，ω1=ω2=1/2

這就是梯形近似的預(yù)估計(jì)-校正法。其他階的龍格-庫塔法可以類推

三階龍格-庫塔法公式：四階龍格-庫塔法公式：例題龍格-庫塔法的近似計(jì)算(4階）解：取步長(zhǎng)h=0.23、有限差分法基礎(chǔ)設(shè)在一個(gè)邊界為S的二維區(qū)域內(nèi)，電位

i,j滿足第一類邊界條件，電位函數(shù)

滿足泊松方程或者拉普拉斯方程。第一類邊界條件：h稱為步長(zhǎng)與節(jié)點(diǎn)(xi，yj)

在x方向直接相鄰的節(jié)點(diǎn)上的電位值表示為將以上兩式相加，并略去h2以上的高階項(xiàng)得與節(jié)點(diǎn)(xi，yj)

在y方向直接相鄰的節(jié)點(diǎn)電位表示為將以上兩式相加，并略去h2以上的高階項(xiàng)得代入拉普拉斯方程代入泊松方程簡(jiǎn)單迭代法對(duì)求解域內(nèi)的節(jié)點(diǎn)賦予迭代初值，可以設(shè)為0。按照如下公式反復(fù)迭代，直到誤差小于給定的迭代閾值。為止。

超松弛迭代法利用超松弛迭代法加快迭代速度，其迭代過程為

為松弛因子。

M為節(jié)點(diǎn)數(shù)。例題

設(shè)有一個(gè)截面為正方形的無限長(zhǎng)接地金屬槽，如圖(a)所示，其導(dǎo)體蓋板電位為10V，并與側(cè)壁絕緣。求其槽內(nèi)電位分布。槽體截面劃分為4×4的16個(gè)正方形網(wǎng)格，共有25個(gè)節(jié)點(diǎn)，其中沿邊界的16個(gè)節(jié)點(diǎn)的電位是已知的，即

所求解的節(jié)點(diǎn)電位為經(jīng)過k=16次的迭代運(yùn)算差分與差商

前向差分后向差分中心差分

4、Poisson方程差分格式

二維Poisson方程如何離散方程的有限差分格式？首先是空間離散，采用什么樣離散？通常采用完全有規(guī)律的分布方式，這樣可使每個(gè)離散點(diǎn)上得到相同形式的差分方程，有效的提高解題速度。對(duì)能填滿平面域的三種規(guī)則網(wǎng)格（正方形，正三角形和正六邊形）的劃分方式，經(jīng)常采用的是矩形網(wǎng)格劃分，Poisson差分格式

一階偏導(dǎo)數(shù)差分格式如何得到高精度的差分格式？可采用待定系數(shù)的分析方法，它的思路：1、3結(jié)點(diǎn)與0結(jié)點(diǎn)在x方向的差分用泰勒公式展開，它們各自占有一定的權(quán)系數(shù)，以截?cái)嗾`差來計(jì)算系數(shù)

Poisson差分格式

一階偏導(dǎo)數(shù)差分格式，忽略h3以上的高次冪的項(xiàng)，并且令項(xiàng)的系數(shù)為零，這樣處理可以保證得到的差分格式誤差為h3量級(jí)。系數(shù)為零的條件

求出具有二階精度精度為一階偏導(dǎo)數(shù)差分格式

Poisson差分格式

二階偏導(dǎo)數(shù)的差分格式，令方程右邊的一階偏導(dǎo)數(shù)的系數(shù)為0，得到系數(shù)間的表達(dá)式代入上式得到精度為O(h3)的二階偏導(dǎo)數(shù)的差分格式

Poisson差分格式

當(dāng)時(shí)，上式可以簡(jiǎn)化為得到Possion方程的五點(diǎn)差分格式5、媒質(zhì)分界面上的差分格式

分界面與網(wǎng)格線重合的情況在0，2，4點(diǎn)的電位值相等兩式中和是假設(shè)“虛”電位，可以利用分界面上場(chǎng)量遵循的邊界條件，削去它們

媒質(zhì)界面的差分格式

由邊界條件，假設(shè)沒有電荷中心差分格式表示把前面關(guān)于和式子代入上式斜媒質(zhì)界面上的差分格式

分界面與網(wǎng)格線呈對(duì)角線的情況在0點(diǎn)上的電位值相同兩式中和是假設(shè)“虛”電位，可以利用分界面上場(chǎng)量遵循的邊界條件，削去它們

不同媒質(zhì)分界面上的差分格式

由邊界條件，假設(shè)界面上沒有電荷對(duì)M、N結(jié)點(diǎn)應(yīng)用線性插值

不同媒質(zhì)分界面上的差分格式

把前面的和代入上式，得網(wǎng)格線呈對(duì)角線的差分格式：對(duì)于矢量磁位的情況：分界面與網(wǎng)格線重合，有拐點(diǎn)（編號(hào)與課本不同）在/2～3

/2，

在0～/2之間變化時(shí)，差分方程和邊界條件不受影響。(1)在0～/2之間變化時(shí)B區(qū)媒質(zhì)滿足拉普拉斯方程：Ab1+Ab2+Ab3+Ab4-4A0=0(2)

=時(shí),退化為分界面與網(wǎng)格線呈對(duì)角線的情況

MN31420B區(qū)A區(qū)得上述（1）為空間只有B區(qū)媒質(zhì)的邊界條件差分格式；（2）為一半空間填充A區(qū)媒質(zhì)的情況下邊界條件差分格式；（B+B/2＋A/2)/2：A/4取兩者的平均，則得到1/4空間填充A區(qū)媒質(zhì)的情況下邊界條件差分格式：也可以利用唯一性定理：（邊界條件不變）（1）兩種情況下，節(jié)點(diǎn)1、4均在B區(qū)，取平均不變；（2）節(jié)點(diǎn)2、3在A、B區(qū)邊界，取平均不變；31420B區(qū)A區(qū)31420B區(qū)A區(qū)+=31420B區(qū)A區(qū)/2根據(jù)得根據(jù)對(duì)稱性，同理兩式相加，除以2媒質(zhì)分界面差分格式（積分法）媒質(zhì)分界面上電位移矢量連續(xù)對(duì)上式在邊界進(jìn)行面積分，利用Gauss定律（注意線積分方向）對(duì)a1,a2,a3,a4四個(gè)邊積分，如a2（只考慮大?。┟劫|(zhì)分界面差分格式（積分法）其他3個(gè)邊類似處理，4個(gè)邊的和整理得6、邊界條件的差分格式第一類邊界條件的差分離散化應(yīng)用多元函數(shù)的泰勒公式，結(jié)點(diǎn)1、3位函數(shù)值通過表示為以h和h1分別與以上兩式相乘且相加，消除一階偏導(dǎo)項(xiàng)，然后截?cái)嗯ch的二次項(xiàng)，便得到關(guān)于結(jié)點(diǎn)0的二階偏導(dǎo)數(shù)的差分格式第一類邊界差分格式同理，在0結(jié)點(diǎn)處關(guān)于y方向的二階偏導(dǎo)的差分格式代入給定的泊松方程，得到通常第一類邊界條件的鄰近點(diǎn)差分格式第三類邊界差分格式第三類邊界條件的差分離散化第一種情況，當(dāng)結(jié)點(diǎn)剛好著落于邊界線L上時(shí)，這還取決于邊界結(jié)點(diǎn)處的外法線與網(wǎng)格線重合，第三類邊界差分格式外法線與網(wǎng)格線不重合情況，邊界結(jié)點(diǎn)上的外向法線方向與水平夾角為ā，其法向?qū)?shù)顯然是在x和y方向的導(dǎo)數(shù)在法向的投影組合，第三類邊界差分格式第二種情況，當(dāng)結(jié)點(diǎn)不落于邊界線L上時(shí)，只需要引入于結(jié)點(diǎn)0相關(guān)的邊界結(jié)點(diǎn)O‘，點(diǎn)的外方向n作為結(jié)點(diǎn)0處的“外方向n”，且近似地認(rèn)為邊界條件中給定的函數(shù)均在O’點(diǎn)上的取值。這樣，此種情況下的第三類邊界條件的離散格式與上式相似，第二類邊界差分格式第二類邊界條件的差分離散化第二類齊次邊界條件為第三類邊界條件的特殊情況，即。我們這里討論最常見的一種情況加一層虛擬邊界上面也是對(duì)稱邊界條件的離散公式第二類邊界差分格式的插值法做法線PQ；OP＝ah，PR＝bh，VP＝ch節(jié)點(diǎn)0與P間滿足差分格式VnPR0B區(qū)A區(qū)chahbhQP點(diǎn)電位由V和R的插值來求出：代入上式，并考慮到第二類邊界差分格式為其中，＝0為第二類邊界的情況，

0為第三類邊界的情況

邊界差分格式可推廣到正三角形六點(diǎn)、六邊形三點(diǎn)式、正方形九點(diǎn)式等情況，見P124直角坐標(biāo)系拉普拉斯方程的邊界的差分格式，見P125拋物型和雙曲型偏微分方程的有限差分法拋物型偏微分方程的有限差分法空間步長(zhǎng)：h=L/J,時(shí)間步長(zhǎng)：τ=T/N1.最簡(jiǎn)單的顯式差分格式令初始條件離散化：左、右端點(diǎn)有關(guān)系：另一方法：泰勒展開的“逐項(xiàng)逼近法”一階展開二階展開利用差分與差商的關(guān)系前項(xiàng)和后項(xiàng)差分為：以上兩式相減得：此即P136的二階泰勒展開的逼近形式。2.最簡(jiǎn)單的隱式差分格式3.六點(diǎn)對(duì)稱格式顯式格式與隱式格式結(jié)合起來表達(dá)了n和n-1層兩層上相鄰6個(gè)結(jié)點(diǎn)的函數(shù)值關(guān)系，故稱為六點(diǎn)差分格式。4.一般線性拋物型微分方程令：①當(dāng)θ=0時(shí)，可以得出顯式格式；②當(dāng)θ≠0時(shí)，得隱式格式。雙曲型偏微分方程的有限差分法，P141顯式格式初速度條件用到了二階差商格式。第一種第二種2.隱式格式第n-1層、n層、n＋1層的中心差商取平均7、差分方程組的求解

綜上所述，對(duì)場(chǎng)域D內(nèi)各結(jié)點(diǎn)（包括所有場(chǎng)域內(nèi)結(jié)點(diǎn)和邊界結(jié)點(diǎn)）逐一列出對(duì)應(yīng)的差分計(jì)算格式，即構(gòu)成以這些離散結(jié)點(diǎn)上的位函數(shù)為待求量的差分方程組（代數(shù)方程組）。求解這些代數(shù)方程組，得到場(chǎng)域中的電位值計(jì)算步驟通常是：離散場(chǎng)域，采用一定的網(wǎng)格剖分方式離散化計(jì)算區(qū)域。離散化場(chǎng)方程，即基于差分原理的應(yīng)用，對(duì)場(chǎng)域內(nèi)場(chǎng)的偏微分方程以及定解條件進(jìn)行差分化處理，得到方程的差分格式。計(jì)算離散解，建立的差分格式（與原定解問題對(duì)立的離散數(shù)學(xué)模型—代數(shù)方程組），選用合適的代數(shù)方程組解法，編寫相應(yīng)的計(jì)算程序，算出待求的結(jié)點(diǎn)上場(chǎng)值。差分法處理過程有限差分法格式特點(diǎn)仔細(xì)分析泊松方程的離散差分方程組，從離散方程式不難看出，該方程組的系數(shù)一般是有規(guī)律的，且方程都很簡(jiǎn)單，每個(gè)方程的項(xiàng)數(shù)不多（待求量最多不超過5項(xiàng)）

各離散結(jié)點(diǎn)上的方程組形式（結(jié)點(diǎn)順序按坐標(biāo)先從y軸增加、再x軸增加（從下到上、從左到右，即先列后行）排列

有限差分法格式特點(diǎn)有限差分法格式特點(diǎn)寫成矩陣方程形式有限差分法格式特點(diǎn)可以看出系數(shù)矩陣由如下特點(diǎn)：系數(shù)矩陣是稀疏矩陣，只有少數(shù)元素不為零。系數(shù)矩陣在一定邊界條件下（邊界與結(jié)點(diǎn)重合且場(chǎng)域邊界類型都一樣），是對(duì)稱正定矩陣。系數(shù)矩陣是方陣，大小為場(chǎng)域中離散結(jié)點(diǎn)的總數(shù)目N=Nx*Ny。，超松弛迭代法計(jì)算上面的線性代數(shù)方程組，通常有兩種方法可以選擇：直接法和迭代法。

松弛算法的基礎(chǔ)為余數(shù)的概念，即余數(shù)逼近零的算法而實(shí)際計(jì)算有一定的誤差超松弛迭代法要計(jì)算出場(chǎng)域內(nèi)的所有結(jié)點(diǎn)上差分格式的解，就要求所有點(diǎn)的差分方程余數(shù)為零。通常情況下，根據(jù)一定的計(jì)算精度要求，各結(jié)點(diǎn)分布的余數(shù)比較均勻且小于要求的精度，就得到滿足精度的近似解。松弛法就是使方程有規(guī)則地朝余數(shù)減少的方向進(jìn)化，直到認(rèn)為得到滿足精度的近似值。為加快松弛算法的計(jì)算速度，科學(xué)家發(fā)展多種松弛法超松弛迭代法求解具有稀疏系數(shù)矩陣的大型差分方程組，其中最優(yōu)的就是超松弛迭代法（SuccessiveOverRelaxation，SOR）。為了說明SOR方法，首先介紹雅可比法和高斯-賽德爾法

雅可比法（Jacobi）就是要使迭代值能精確的滿足前一次各點(diǎn)的電位值所能表示的差分方程超松弛迭代法高斯-賽德爾法是雅可比法的改進(jìn)方法，主要針對(duì)減少內(nèi)存消耗，只需存儲(chǔ)一組完整的數(shù)組。它采取的措施是對(duì)每一次迭代盡量采用最新計(jì)算的值來替換上一次迭代的舊值。結(jié)果收斂速度比雅可比法快一倍。超松弛迭代法逐次超松弛法是對(duì)高斯-賽德爾法的改進(jìn)，該方法的核心是借助于一收斂因子ω作用到高斯-賽德爾迭代公式。當(dāng)時(shí)w＝1，就回到高斯-賽德爾法。當(dāng)w>2時(shí)，迭代過程變得及其不穩(wěn)定。只有1<w<2，才能提高收斂速度。

超松弛迭代法正方形第一類邊界條件時(shí)長(zhǎng)方形第一類邊界條件時(shí)場(chǎng)強(qiáng)與電、磁積分量的計(jì)算通過上述差分方程組的求解，在獲得場(chǎng)域內(nèi)各結(jié)點(diǎn)上待求位函數(shù)后，往往還需求場(chǎng)中的場(chǎng)強(qiáng)分布，以及其他有關(guān)的積分特性（如磁通量和磁導(dǎo)、電導(dǎo)、電容等磁路及電路參數(shù)等）。

場(chǎng)強(qiáng)與電、磁積分量的計(jì)算無論是靜電場(chǎng)、恒定電流場(chǎng)或恒定磁場(chǎng)，其通量可一般地表示為所分析的靜電場(chǎng)中的電容C、恒定電流場(chǎng)中的電導(dǎo)G或恒定磁場(chǎng)中的磁導(dǎo)等電路或磁路參數(shù)P就可按下式計(jì)算收斂因子.與收斂性的關(guān)系一個(gè)穩(wěn)定的收斂過程，降低任一結(jié)點(diǎn)最大誤差到預(yù)定的誤差

所需要的迭代次數(shù)N由下式?jīng)Q定:式中:F為漸近收斂速度，是邊界形狀、邊界條件、結(jié)點(diǎn)數(shù)、差分方程形式、收斂因子的函數(shù)。有限差分法的收斂性和穩(wěn)定性下面僅就單步法進(jìn)行探討。尤拉折線法、龍格一庫塔法、泰勒展開法等:單步法的一般形式可以表示為單步法的截?cái)嗾`差為一般，單步法的截?cái)嗾`差可以用泰勒展開來估算稱為k階方法在y滿足李普西茲(Iipschits)條件;計(jì)算誤差可用