遼工大有限差分法實(shí)驗(yàn)報(bào)告

1、電磁場(chǎng)與電磁波項(xiàng)目訓(xùn)練報(bào)告求解金屬槽的電位分布班 級(jí)： 姓 名： 學(xué) 號(hào)： 指導(dǎo)教師： 成 績(jī)： 電子與信息工程學(xué)院信息與通信工程系求解金屬槽的電位分布1.實(shí)驗(yàn)原理利用有限差分法和matlab軟件解決電位在金屬槽中的分布。有限差分法基本思想是把連續(xù)的定解區(qū)域用有限個(gè)離散點(diǎn)構(gòu)成的網(wǎng)格來(lái)代替,這些離散點(diǎn)稱(chēng)作網(wǎng)格的節(jié)點(diǎn)；把連續(xù)定解區(qū)域上的連續(xù)變量的函數(shù)用在網(wǎng)格上定義的離散變量函數(shù)來(lái)近似；把原方程和定解條件中的微商用差商來(lái)近似,積分用積分和來(lái)近似,于是原微分方程和定解條件就近似地代之以代數(shù)方程組,即有限差分方程組 ,解此方程組就可以得到原問(wèn)題在離散點(diǎn)上的近似解.然后再利用插值方法便可以從離散解得到定解

2、問(wèn)題在整個(gè)區(qū)域上的近似解.在采用數(shù)值計(jì)算方法求解偏微分方程時(shí),若將每一處導(dǎo)數(shù)由有限差分近似公式替代,從而把求解偏微分方程的問(wèn)題轉(zhuǎn)換成求解代數(shù)方程的問(wèn)題。2.有限差分法方程的定解問(wèn)題就是在滿(mǎn)足某些定解條件下求微分方程的解。在空間區(qū)域的邊界上要滿(mǎn)足的定解條件稱(chēng)為邊值條件。如果問(wèn)題與時(shí)間有關(guān)，在初始時(shí)刻所要滿(mǎn)足的定解條件，稱(chēng)為初值條件。不含時(shí)間而只帶邊值條件的定解問(wèn)題，稱(chēng)為邊值問(wèn)題。與時(shí)間有關(guān)而只帶初值條件的定解問(wèn)題，稱(chēng)為初值問(wèn)題。同時(shí)帶有兩種定解條件的問(wèn)題，稱(chēng)為初值邊值混合問(wèn)題。定解問(wèn)題往往不具有解析解，或者其解析解不易計(jì)算。所以要采用可行的數(shù)值解法。有限差分方法就是一種數(shù)值解法，它的基本思想是先

3、把問(wèn)題的定義域進(jìn)行網(wǎng)格剖分，然后在網(wǎng)格點(diǎn)上，按適當(dāng)?shù)臄?shù)值微分公式把定解問(wèn)題中的微商換成差商，從而把原問(wèn)題離散化為差分格式，進(jìn)而求出數(shù)值解。此外，還要研究差分格式的解的存在性和唯一性、解的求法、解法的數(shù)值穩(wěn)定性、差分格式的解與原定解問(wèn)題的真解的誤差估計(jì)、差分格式的解當(dāng)網(wǎng)格大小趨于零時(shí)是否趨于真解（即收斂性），等等。有限差分方法具有簡(jiǎn)單、靈活以及通用性強(qiáng)等特點(diǎn)，容易在計(jì)算機(jī)上實(shí)現(xiàn)。2.1有限差分法原理 圖1-1 有限差分法的網(wǎng)格劃分導(dǎo)體槽中靜電場(chǎng)的邊值問(wèn)題的拉普拉斯方程為： (1-1)為簡(jiǎn)單起見(jiàn)，將場(chǎng)域分成足夠小的正方形網(wǎng)格，網(wǎng)格線(xiàn)之間的距離為h，。節(jié)點(diǎn)0、1、2、3、4上的電位分別用、和表示。點(diǎn)

4、1、點(diǎn)3在x0處可微，沿x方向在x0處的泰勒級(jí)數(shù)展開(kāi)式為 (1-2) (1-3)點(diǎn)2、點(diǎn)4在y0處可微，沿y方向在y0處的泰勒級(jí)數(shù)展開(kāi)式為 (1-4) (1-5)忽略高次項(xiàng) (1-6)稍作變化得到拉普拉斯方程的五點(diǎn)差分格式： (1-7) 可通過(guò)迭代法求解以上差分方程。2.2有限差分法步驟高斯賽德?tīng)柕▓D2-1網(wǎng)絡(luò)下標(biāo)標(biāo)示 (1-8)進(jìn)行迭代時(shí)可寫(xiě)為 (1-9),為行數(shù)，,為列數(shù)，為迭代次數(shù)，為前次迭代的結(jié)果，為當(dāng)次迭代的結(jié)果，由于迭代從第一行、第一列開(kāi)始，（）、（）點(diǎn)的迭代較（）點(diǎn)進(jìn)行得早，顧可使用當(dāng)次迭代的結(jié)果。直到所有的點(diǎn)電位滿(mǎn)足（為所設(shè)定精度）時(shí)迭代停止。3.問(wèn)題描述設(shè)有一個(gè)長(zhǎng)直接地金

5、屬矩形槽，如圖2-1所示，其側(cè)壁與底面點(diǎn)位均為零，頂蓋電位為100V（相對(duì)值），求槽內(nèi)點(diǎn)位分布。圖3-1 金屬槽4.程序設(shè)計(jì)4.1全場(chǎng)域問(wèn)題對(duì)于問(wèn)題（1）（2）（3）而言，以步距的正方形網(wǎng)格離散化場(chǎng)域。每個(gè)網(wǎng)格對(duì)應(yīng)于矩陣中的單個(gè)元素。由此通過(guò)矩陣中的值的計(jì)算并指定相鄰兩次的迭代值誤差不超過(guò)，應(yīng)用matlab中的矩陣操作。利用ones(x,y)建立一個(gè)且每個(gè)元素初值都為為1的矩陣A。再對(duì)矩陣進(jìn)行狄利克雷邊界初始化，并且設(shè)置矩陣的左右邊界為0，上下邊界分別為100和0。在保證精度的情況下以拉普拉斯差分式進(jìn)行下一級(jí)的數(shù)值計(jì)算。最終得到一個(gè)滿(mǎn)足迭代要求的矩陣A。具體程序?qū)崿F(xiàn)見(jiàn)附錄A。4.2程序流程圖圖

6、4-1 程序設(shè)計(jì)流程圖通過(guò)函數(shù)contour(A)可以繪制出等電位分布圖，這樣可以觀察出與理論情況的分布是否相同，分布圖如圖2-3所示：圖4-2 金屬槽等電位圖4.3半場(chǎng)域問(wèn)題對(duì)于問(wèn)題（4），在程序設(shè)計(jì)中利用ones(x,y)建立一個(gè)且每個(gè)元素初值都為為1的矩陣A。再對(duì)矩陣進(jìn)行狄利克雷邊界初始化，并且設(shè)置矩陣的左邊界為0，左邊界為初值為由式子100:-5:0產(chǎn)生的一個(gè)單位矩陣向量。上下邊界分別為100和0。在保證精度的情況下以拉普拉斯差分式進(jìn)行下一級(jí)的數(shù)值計(jì)算。最終得到一個(gè)滿(mǎn)足迭代要求的矩陣A。具體程序?qū)崿F(xiàn)見(jiàn)附錄B。運(yùn)行附錄B中的程序就可以得到半場(chǎng)域的數(shù)值解。通過(guò)函數(shù)contour(A)可以繪

7、制出等電位分布圖，分布圖如下：圖4-3 半場(chǎng)域金屬槽等電位圖4.4中心點(diǎn)處的討論對(duì)于問(wèn)題（7），取中心點(diǎn)P()。求解其電位的解析解，因?yàn)樵擖c(diǎn)位置處于高度對(duì)稱(chēng)位置，所以該點(diǎn)所在的對(duì)稱(chēng)線(xiàn)上可視為勻強(qiáng)電場(chǎng)。故有：E · L = (4-1)得： (4-2)所以該點(diǎn)電位。由附錄A中的程序運(yùn)行得出該點(diǎn)的數(shù)值解為49.516V。由此可以得出誤差范圍為：0.4840V，相對(duì)誤差為：0.9700%。5.收斂因子作用5.1超松弛迭代法為了加快收斂速度，常采用超松弛迭代法。計(jì)算時(shí)，將某點(diǎn)的新老電位值之差乘以一個(gè)因子以后，再加到該點(diǎn)的新老電位值上，作為這一點(diǎn)的新電位值。超松弛迭代法的表達(dá)式： (5-1)式中

8、稱(chēng)為松弛因子，其值介于1和2之間。其中最佳收斂因子: (5-2)其中為每邊的節(jié)點(diǎn)數(shù)減去1。 其中M、N分別是x、y兩個(gè)方向的內(nèi)節(jié)點(diǎn)數(shù)。對(duì)于本項(xiàng)目中，M=41、N=21，計(jì)算得：表4-1 最佳收斂因子迭代次數(shù)1.7516迭代次數(shù)89取若干個(gè)收斂因子求得的迭代次數(shù)得表4-2和表4-3中的數(shù)據(jù)。表4-2 收斂因子迭代次數(shù)和中心數(shù)值1.01.11.21.31.41.5迭代次數(shù)740615508416335262中心點(diǎn)電位值49.516749.516849.516949.517049.517149.5171表4-3 收斂因子迭代次數(shù)和中心數(shù)值1.61.71.741.781.81.9迭代次數(shù)1961331

9、067380161中心點(diǎn)電位值49.517249.517249.517349.517349.517349.5173其中：中心點(diǎn)的解析解為50V。由收斂因子映射到中心電位值可以得到一個(gè)利用Excel繪制的走勢(shì)圖。如圖：圖4-1 中心點(diǎn)電勢(shì)趨勢(shì)圖5.2關(guān)于收斂因子的結(jié)論5.2.1從迭代次數(shù)角度分析從迭代次數(shù)上看，收斂因子存在一個(gè)最佳值。由表可以看出，當(dāng)?shù)蜃訌?開(kāi)始趨近于最佳迭代因子時(shí)，收斂次數(shù)減少，進(jìn)而使收斂速度較簡(jiǎn)單迭代法加快。當(dāng)收斂因子遠(yuǎn)離最佳收斂因子時(shí)，收斂次數(shù)又開(kāi)始減少，進(jìn)而使收斂速度較簡(jiǎn)單迭代法變慢。即收斂因子影響收斂速度。5.2.2從迭代次數(shù)角度分析從數(shù)值解的角度來(lái)說(shuō)，由圖3-7可

10、知，隨著加速收斂因子的增加，所得到的數(shù)值解就越接近于也隨之增加。到達(dá)一定的數(shù)值后，呈現(xiàn)出一種穩(wěn)定的狀態(tài)。由此可知，收斂因子可以的增加會(huì)使數(shù)值解接近于解析解，而又不會(huì)等于解析解。參考文獻(xiàn)1 司守奎，孫璽箐.數(shù)學(xué)建模算法與應(yīng)用.國(guó)防工業(yè)出版社，2015.2:4114242 王家禮，朱滿(mǎn)座，路宏敏.電磁場(chǎng)與電磁波.西安電子科技大學(xué)出版社，2009.8:118122附錄A：簡(jiǎn)單差分法求解電位分布%將待求區(qū)域化成20*40個(gè)邊長(zhǎng)為1的正方體% 頂板100V，左右底為0Vx=21;y=41;%網(wǎng)格節(jié)點(diǎn)數(shù) x行，y列A=ones(x,y);%設(shè)置21行，41列的矩陣A(1,:)=ones(1,y)*100;

11、% 設(shè)置上基板的值A(chǔ)(x,:)=zeros(1,y);% 設(shè)置下基板的值A(chǔ)(:,1)=zeros(x,1);% 設(shè)置左基板的值A(chǔ)(:,y)=zeros(x,1);% 設(shè)置右基板的值disp(A);%命令行輸出矩陣AA1=A;%初始化一級(jí)近似值A(chǔ)1max=1;%初始化最大絕對(duì)誤差值，用于進(jìn)入while循環(huán)k=0;%迭代次數(shù)while(max>1e-5)%由A迭代，算出·A1，迭代精度為1e-6 k=k+1;%計(jì)算迭代次數(shù) max=0;%誤差回歸最小 temp=0;%單次計(jì)算的前后兩次迭代中的同元素的誤差初始值 for i=2:1:x-1 %從第2行到底20行 for j=2:1:

12、y-1 %從第2列到底40列 A1(i,j)=( A(i,j+1)+A(i,j-1)+A(i+1,j)+A(i-1,j) )/4;%拉普拉斯方程差分式 disp('A1(i,j)='),disp(A1(i,j);%log輸出 temp=abs(A1(i,j)-A(i,j);%相鄰兩次迭代解之間的誤差 disp('temp='),disp(temp);%log輸出 %控制精度的最大值，得到誤差計(jì)算中的最大值 if temp>max max=temp; end end end A=A1;%賦值給下次用enddisp('輸出A1='),disp(

13、A);disp('輸出迭代次數(shù)k='),disp(k);附錄B：收斂因子作用程序設(shè)計(jì)%比較最佳收斂因子和收斂因子的差異x=21;y=41;%網(wǎng)格節(jié)點(diǎn)數(shù) x行，y列A=ones(x,y);%設(shè)置21行，41列的矩陣A(1,:)=ones(1,y)*100;% 設(shè)置上基板的值A(chǔ)(x,:)=zeros(1,y);% 設(shè)置下基板的值A(chǔ)(:,1)=zeros(x,1);% 設(shè)置左基板的值A(chǔ)(:,y)=zeros(x,1);% 設(shè)置右基板的值disp(A);%命令行輸出矩陣A瞅瞅A1=A;%初始化一級(jí)近似值A(chǔ)1max=1;%初始化最大絕對(duì)誤差值，用于進(jìn)入while循環(huán)k=0;%迭代次數(shù)a0

14、=2-pi*sqrt(2/(412)+(2/(212);%最佳收斂因子%-注意事項(xiàng)-%修改下面a的值就可以任意設(shè)定收斂因子%取值需要滿(mǎn)足：1<=a<2%-注意事項(xiàng)-a=a0 ;% 收斂因子while(max>1e-5)%由A迭代，算出·A1，迭代精度為1e-6 k=k+1;%計(jì)算迭代次數(shù) max=0;%誤差回歸最小 temp=0;%單次計(jì)算的前后兩次迭代中的同元素的誤差初始值 for i=2:1:x-1 %從第2行到底20行 for j=2:1:y-1 %從第2列到底40列 A1(i,j)=A(i,j)+( A(i,j+1)+A1(i,j-1)+A(i+1,j)+A

15、1(i-1,j)-4*A(i,j) )*a/4;%超松弛迭代法表達(dá)式 disp('A1(i,j)='),disp(A1(i,j);%log輸出 temp=abs(A1(i,j)-A(i,j);%相鄰兩次迭代解之間的誤差 disp('temp='),disp(temp);%log輸出 %控制精度的最大值，得到誤差計(jì)算中的最大值 if temp>max max=temp; end end end A=A1;%賦值給下次用enddisp('輸出A1='),disp(A1);disp('輸出迭代次數(shù)k='),disp(k);附錄C：

16、半場(chǎng)域程序設(shè)計(jì)%半場(chǎng)域的計(jì)算x=21;y=21;%網(wǎng)格節(jié)點(diǎn)數(shù) x行，y列A=ones(x,y);%設(shè)置21行，21列的矩陣A(1,:)=ones(1,y)*100;% 設(shè)置上基板的值A(chǔ)(x,:)=zeros(1,y);% 設(shè)置下基板的值A(chǔ)(:,1)=zeros(x,1);% 設(shè)置左基板的值right=100:-5:0;%建立一個(gè)100開(kāi)頭的行向量A(:,y)=right;% 設(shè)置右基板的值disp(A);%命令行輸出矩陣AA1=A;%初始化一級(jí)近似值A(chǔ)1max=1;%初始化最大絕對(duì)誤差值，用于進(jìn)入while循環(huán)k=0;%迭代次數(shù)while(max>1e-5)%由A迭代，算出·A1，迭代精度為1e-6 k=k+1;%計(jì)算迭代次數(shù) max=0;%誤差回歸最小 temp=0;%單次計(jì)算的前后兩次迭代中的同元素的誤差初始值 for i=2:1:x-1 %從第2行到底21行 for j=2:1:y-1 %從第2列到底21列 A1(i,j)=( A(i,j+1)+A(i,j-1)+A(i+1,j)+A(i-1,j) )/4;%拉普拉斯方程差分式 d