有限差分法

班级：通信13-4 姓名：
学号：
指导教师：**
成绩：
电子与信息工程学院
信息与通信工程系
求解金属槽的电位分布
1.实验原理
利用有限差分法和matlab软件解决电位在金属槽中的分布。

有限差分法基本思想是把连续的定解区域用有限个离散点构成的网格来代替,这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似,积分用积分和来近似,于是原微分方程和定解条件就近似地代之以代数方程组,即有限差分方程组,解此方程组就可以得到原问题在离散点上的近似解.然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解.在采用数值计算方法求解偏微分方程时,若将每一处导数由有限差分近似公式替代,从而把求解偏微分方程的问题转换成求解代数方程的问题。

2.有限差分法
方程的定解问题就是在满足某些定解条件下求微分方程的解。

在空间区域的边界上要满足的定解条件称为边值条件。

如果问题与时间有关，在初始时刻所要满足的定解条件，称为初值条件。

不含时间而只带边值条件的定解问题，称为边值问题。

与时间有关而只带初值条件的定解问题，称为初值问题。

同时带有两种定解条件的问题，称为初值边值混合问题。

定解问题往往不具有解析解，或者其解析解不易计算。

所以要采用可行的数值解法。

有限差分方法就是一种数值解法，它的基本思想是先把问题的定义域进行网格剖分，然后在网格点上，按适当的数值微分公式把定解问题中的微商换成差商，从而把原问题离散化为差分格式，进而求出数值解。

此外，还要研究差分格式的解的存在性和唯一性、解的求法、解法的数值稳定性、差分格式的解与原定解问题的真解的误差估计、差分格式的解当网格大小趋于零时是否趋于真解（即收敛性），等等。

有限差分方法具有简单、灵活以及通用性强等特点，容易在计算机上实现。

2.1有限差分法原理
图1-1 有限差分法的网格划分
导体槽中静电场的边值问题的拉普拉斯方程为：
22
22
0x y ϕϕ∂∂+=∂∂ (1-1) 为简单起见，将场域分成足够小的正方形网格，网格线之间的距离为h ，0h →。

节点0、1、2、3、4上的电位分别用0ϕ、1ϕ、2ϕ、3ϕ和4ϕ表示。

点1、点3在x 0处可微，沿x 方向在x 0处的泰勒级数展开式为
23
23100002311()()()()().2!3!h h h x x x ϕϕϕϕϕ∂∂∂=-+-+-+∂∂∂ (1-2)
2323300002311()()()()().2!3!h h h x x x
ϕϕϕϕϕ∂∂∂=-+-+-+∂∂∂ (1-3)
点2、点4在y 0处可微，沿y 方向在y 0处的泰勒级数展开式为
2323200002311()()().
2!3!h h h y y y
ϕϕϕ
ϕϕ∂∂∂=++++∂∂∂ (1-4)
2323400002311()()()()().2!3!h h h y y y
ϕϕϕϕϕ∂∂∂=-+-+-+∂∂∂ (1-5)
忽略高次项
222
12340000224()()()4h x
y ϕϕϕϕϕϕϕϕ⎡⎤
∂∂+++=++=⎢⎥∂∂⎣⎦ (1-6)
稍作变化得到拉普拉斯方程的五点差分格式：
1234
04
ϕϕϕϕϕ+++=
(1-7)
可通过迭代法求解以上差分方程。

2.2有限差分法步骤
高斯—赛德尔迭代法
图2-1网络下标标示
1,,1,11,,4
i j i j i j i j
i j ϕϕϕϕϕ--+++++=
(1-8)
进行迭代时可写为 ,1
,111,,11,14
i j i j k k k k i j i j i j k ϕϕ
ϕϕϕ-++-++++++=
(1-9)
1,2.......i M =,为行数，1,2.......j N =,为列数，k 为迭代次数，k ϕ为前次迭代的结果，1k ϕ+为
当次迭代的结果，由于迭代从第一行、第一列开始，（1,i j -）、（,1i j -）点的迭代较（,i j ）点进行得早，顾可使用当次迭代的结果。

直到所有的点电位满足,,1i j i j k k ϕϕε+-<（ε为所设定精度）时迭代停止。

3.问题描述
设有一个长直接地金属矩形槽，如图2-1所示，其侧壁与底面点位均为零，顶盖电位为100V （相对值），求槽内点位分布。

图3-1 金属槽
4.程序设计
4.1全场域问题
对于问题（1）（2）（3）而言，以步距的正方形网格离散化场域。

每个网格对应于矩阵中的单个元素。

由此通过矩阵中的值的计算并指定相邻两次的迭代值误差不超过，应用matlab中的矩阵操作。

利用ones(x,y)建立一个且每个元素初值都为为1的矩阵A。

再对矩阵进行狄利克雷边界初始化，并且设置矩阵的左右边界为0，上下边界分别为100和0。

在保证精度的情况下以拉普拉斯差分式进行下一级的数值计算。

最终得到一个满足迭代要求的矩阵A。

具体程序实现见附录A。

4.2程序流程图
图4-1 程序设计流程图
通过函数contour(A)可以绘制出等电位分布图，这样可以观察出与理论情况的分布是否相同，分布图如图2-3所示：
图4-2 金属槽等电位图
4.3半场域问题
对于问题（4），在程序设计中利用ones(x,y)建立一个且每个元素初值都为为1的矩阵A。

再对矩阵进行狄利克雷边界初始化，并且设置矩阵的左边界为0，左边界为初值为由式子100:-5:0产生的一个单位矩阵向量。

上下边界分别为100和0。

在保证精度的情况下以拉普拉斯差分式进行下一级的数值计算。

最终得到一个满足迭代要求的矩阵A。

具体
程序实现见附录B 。

运行附录B 中的程序就可以得到半场域的数值解。

通过函数contour(A)可以绘制出等电位分布图，分布图如下：
图4-3 半场域金属槽等电位图
4.4中心点处的讨论
对于问题（7），取中心点P()。

求解其电位的解析解，因为该点位置处于高度对称位置，所以该点所在的对称线上可视为匀强电场。

故有：
E · L = φ (4-1)
得：
φ/2 = L/2 · E (4-2) 所以该点电位。

由附录A 中的程序运行得出该点的数值解为49.516V 。

由此可以得出误差范围为：0.4840V ，相对误差为：0.9700%。

5.收敛因子作用
5.1超松弛迭代法
为了加快收敛速度，常采用超松弛迭代法。

计算时，将某点的新老电位值之差乘以一个因子以后，再加到该点的新老电位值上，作为这一点的新电位值。

超松弛迭代法的表达式：
,1,,,111,,11,1(4)
4
i j i j i j
i j
k k k k k
i j i j i j k k ϕϕϕϕϕϕϕα
-++-++++++-=+ (5-1)
式中α称为松弛因子，其值介于1和2之间。

其中最佳收敛因子:
2
1sin()
o p απ=
+ (5-2)
其中p 为每边的节点数减去1。

其中M 、N 分别是x 、y 两个方向的内节点数。

对于本项目中，M=41、N=21，计算得：
表4-1 最佳收敛因子迭代次数
1.7516 迭代次数
89
取若干个收敛因子求得的迭代次数得表4-2和表4-3中的数据。

表4-2 收敛因子迭代次数和中心数值
1.0 1.1 1.2 1.3 1.4 1.5 迭代次数 740 615 508 416 335 262 中心点电位值
49.5167
49.5168
49.5169
49.5170
49.5171
49.5171
表4-3 收敛因子迭代次数和中心数值
1.6 1.7 1.74 1.78 1.8 1.9 迭代次数 196 133 106 73 80 161 中心点电位值
49.5172
49.5172
49.5173
49.5173
49.5173
49.5173
其中：中心点的解析解为50V 。

由收敛因子映射到中心电位值可以得到一个利用Excel 绘制的走势图。

如图：
图4-1 中心点电势趋势图
5.2关于收敛因子的结论
5.2.1从迭代次数角度分析
从迭代次数上看，收敛因子存在一个最佳值。

由表可以看出，当迭代因子从1开始趋近于最佳迭代因子时，收敛次数减少，进而使收敛速度较简单迭代法加快。

当收敛因子远离最佳收敛因子时，收敛次数又开始减少，进而使收敛速度较简单迭代法变慢。

即收敛因子影响收敛速度。

5.2.2从迭代次数角度分析
从数值解的角度来说，由图3-7可知，随着加速收敛因子的增加，所得到的数值解就越接近于也随之增加。

到达一定的数值后，呈现出一种稳定的状态。

由此可知，收敛因子可以的增加会使数值解接近于解析解，而又不会等于解析解。

参考文献
[1] 司守奎，孙玺箐.数学建模算法与应用.国防工业出版社，2015.2:411—424
[2] 王家礼，朱满座，路宏敏.电磁场与电磁波.西安电子科技大学出版社，2009.8:118—122
附录A：简单差分法求解电位分布
%将待求区域化成20*40个边长为1的正方体
% 顶板100V，左右底为0V
x=21;y=41;%网格节点数x行，y列
A=ones(x,y);%设置21行，41列的矩阵
A(1,:)=ones(1,y)*100;% 设置上基板的值
A(x,:)=zeros(1,y);% 设置下基板的值
A(:,1)=zeros(x,1);% 设置左基板的值
A(:,y)=zeros(x,1);% 设置右基板的值
disp(A);%命令行输出矩阵A瞅瞅
A1=A;%初始化一级近似值A1
max=1;%初始化最大绝对误差值，用于进入while循环
k=0;%迭代次数
while(max>1e-5)%由A迭代，算出·A1，迭代精度为1e-6
k=k+1;%计算迭代次数
max=0;%误差回归最小
temp=0;%单次计算的前后两次迭代中的同元素的误差初始值
for i=2:1:x-1 %从第2行到底20行
for j=2:1: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输出
%控制精度的最大值，得到误差计算中的最大值
if temp>max
max=temp;
end
end
end
A=A1;%赋值给下次用
end
disp('输出A1===='),disp(A);
disp('输出迭代次数k===='),disp(k);
附录B：收敛因子作用程序设计
%比较最佳收敛因子和收敛因子的差异
x=21;y=41;%网格节点数x行，y列
A=ones(x,y);%设置21行，41列的矩阵
A(1,:)=ones(1,y)*100;% 设置上基板的值
A(x,:)=zeros(1,y);% 设置下基板的值
A(:,1)=zeros(x,1);% 设置左基板的值
A(:,y)=zeros(x,1);% 设置右基板的值
disp(A);%命令行输出矩阵A瞅瞅
A1=A;%初始化一级近似值A1
max=1;%初始化最大绝对误差值，用于进入while循环
k=0;%迭代次数
a0=2-pi*sqrt((2/(41^2))+(2/(21^2)));%最佳收敛因子
%------------------------------注意事项----------------------------
%修改下面a的值就可以任意设定收敛因子
%取值需要满足：1<=a<2
%------------------------------注意事项----------------------------
a=a0 ;% 收敛因子
while(max>1e-5)%由A迭代，算出·A1，迭代精度为1e-6
k=k+1;%计算迭代次数
max=0;%误差回归最小
temp=0;%单次计算的前后两次迭代中的同元素的误差初始值
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)+A1(i-1,j)-4*A(i,j) )*a/4;%超松弛迭代法表达式disp('A1(i,j)='),disp(A1(i,j));%log输出
temp=abs(A1(i,j)-A(i,j));%相邻两次迭代解之间的误差
disp('temp='),disp(temp);%log输出
%控制精度的最大值，得到误差计算中的最大值
if temp>max
max=temp;
end
end
end
A=A1;%赋值给下次用
end
disp('输出A1===='),disp(A1);
disp('输出迭代次数k===='),disp(k);
附录C：半场域程序设计
%半场域的计算
x=21;y=21;%网格节点数x行，y列
A=ones(x,y);%设置21行，21列的矩阵
A(1,:)=ones(1,y)*100;% 设置上基板的值
A(x,:)=zeros(1,y);% 设置下基板的值
A(:,1)=zeros(x,1);% 设置左基板的值
right=100:-5:0;%建立一个100开头的行向量
A(:,y)=right;% 设置右基板的值
disp(A);%命令行输出矩阵A瞅瞅
A1=A;%初始化一级近似值A1
max=1;%初始化最大绝对误差值，用于进入while循环
k=0;%迭代次数
while(max>1e-5)%由A迭代，算出·A1，迭代精度为1e-6
k=k+1;%计算迭代次数
max=0;%误差回归最小
temp=0;%单次计算的前后两次迭代中的同元素的误差初始值
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;%拉普拉斯方程差分式
disp('A1(i,j)='),disp(A1(i,j));%log输出
temp=abs(A1(i,j)-A(i,j));%相邻两次迭代解之间的误差
disp('temp='),disp(temp);%log输出
%控制精度的最大值，得到误差计算中的最大值
if temp>max
max=temp;
end
end
end
A=A1;%赋值给下次用
end
disp('输出A1===='),disp(A1);
disp('输出迭代次数k===='),disp(k);。