静电场边值问题仿真实验-实验指导书

A(k,k+1)=1.0/dx^2;
A(k,k-(m-1))=1.0/dy^2;
A(k,k+(m-1))=1.0/dy^2;
elseif i==1
A(k,k)=-2.0*(1.0/dx^2+1.0/dy^2);
A(k,k+1)=1.0/dx^2;
A(k,k-(m-1))=1.0/dy^2;
A(k,k+(m-1))=1.0/dy^2;
x=0
x=a
Δx
待求区域网格划分
则当 Δx 、 Δy 足够小时，在计算区域内部任一节点 (xi , y j ) 处，有
ϕ(xi + Δx, y j ) − 2ϕ(xi , y j ) + ϕ(xi − Δx, y j ) Δx 2
+ ϕ(xi , y j + Δy) − 2ϕ(xi , y j ) + ϕ(xi , y j − Δy) = 0 Δy 2
实验、二维静电场边值问题的数值仿真实验
对于静电场边值问题，某些情况下用解析法求解很困难，甚至不可能。在实际 求解过程中，可以将待求解电磁场问题化为离散的数值问题，即求解待求量在一些 离散点上的数值。本实验根据二维静电场所满足的微分方程和边界条件，利用有限 差分法来求解无源区域内一系列离散点上的电位，从而观测在给定边界条件下无源 区域静电场的电位分布。
即确定出在边界节点处的值。
第一类边界条件ϕ(x, y) = g(x, y) ：
对于边界上的任一节点 (xi , y j ) ，由于数 g(x, y) 在该点的值 g(xi , y j ) 为已知，可
取
ϕ (xi , y j ) = g(xi , y j )
(8)
这样就给出在边界每一节点上的离散边界条件。
或 Surf 函数）。
2、其它条件不变，在
y=b
的边界上令
f
(x)
=
10
π sin(
x) ，计算并绘制出区域内
a
的电位分布。观察两种情况下电位分布的差别，理解边界条件对场分布的决定作用。
3、对上述两种情况，分别以列表的形式给出 x = 10Δx 和 y = 5Δy 处各节点的电
位与分离变量法（无穷级数截断到第 100 项）所得结果的比较。
% f(i*dx)是上边界给定的已知电位
end
5
end
生成 A 和 b 后，待求向量 V 可以由 V = A −1b 求得。
三、实验内容
1、自行给定区域尺寸 a、b 和，取 Δx = a 20 ，Δy = b 10 ，在 y=b 的边界上令 f(x)=10，
计算各内部节点处的电位值，并绘制出区域内的电位分布（可用 Matlab 中的 Mesh
Δx 2
Δy 2
1
近似地代替。 这里希望根据边界条件来确定矩形区域内的 u(x, y) ，首先，作平行于坐标轴的
两组直线 xi = iΔx y j = jΔy
这两族直线将区域分割成有限数目的网格，网格的边长 Δx 、Δy 称为步长，点 (xi , y j ) 称为网格的节点，如图所示。
y=b
Δy
y=0
（2）
u(x − Δx) = u(x) − u′(x)Δx + 1 u′′(x)(Δx)2 + " 2
（2）+（3）并忽略高次项得：
（3）
u′′(x) ≈ u(x + Δx) − 2u(x) + u(x − Δx) Δx 2
则二维拉普拉斯方程可以用方程
（5）
ϕ(x + Δx, y) − 2ϕ(x, y) + ϕ(x − Δx, y) + ϕ(x, y + Δy) − 2ϕ(x, y) + ϕ(x, y − Δy) = 0 （6）
一、实验目的
1、学习用有限差分法求解静电场的方法。 2、观察边界条件对场的影响。
二、实验原理
1、方程的离散 有限差分法对于解决任何偏微分方程都是一种可行的方法，因为所有的电磁场
问题都可以表示成标量或矢量偏微分方程，所以可以利用有限差分法解决各种介质 中电磁场的空间分布。有限差分法把要解决的区域划分为有限个离散点，并将偏微 分方程用一组差分方程代替。因此，这种方法是近似的。但是，如果我们把离散点 取得足够密，就能够把误差减小到可以接受的程度。
b=zeros((N-1)*(M-1));
% 事先把已知列向量各元素置0，在重新赋值前保持为0
for j=1:N-1
for i=1:M-1
k=(j-1)*(M-1)+i;
if j~=1 & j~=N-1
if i~=1 & i~=M-1
A(k,k-1)=1.0/dx^2;
A(k,k)=-2.0*(1.0/dx^2+1.0/dy^2);
(10e)
经整理，（10）中的各式分别可写为：
[ ] [ ] φk+1 − 2φk + φk−1 Δx2 + φk+(m−1) − 2φk + φk−(m−1) Δy2 = 0
(j≠1、N-1，i≠1、M-1)
[ ] [ ] φk+1 − 2φk Δx2 + φk+(m−1) − 2φk + φk−(m−1) Δy 2 = 0
(j≠1、N-1，i=1)
[ ] [ ] − 2φk + φk−1 Δx2 + φk+(m−1) − 2φk + φk−(m−1) Δy2 = 0
(j≠1、N-1，i=M-1)
[ ] [ ] φk+1 − 2φk + φk−1 Δx2 + φk+(m−1) − 2φk Δy 2 = 0
(j=1)
4
[ ] [ ] φk+1 − 2φk + φk−1 Δx2 + − 2φk + φk−(m−1) Δy2 = − f (iΔx) Δy2
(10c)
当 j=1 时，有
[ ] [ ] φk +1 − 2φk + φk −1 Δx2 + φk +(m−1) − 2φk + 0 Δy2 = 0
(10d)
当 j=N-1 时，有
[ ] [ ] φk+1 − 2φk + φk−1 Δx2 + f (iΔx) − 2φk + φk−(m−1) Δy2 = 0
ϕ(xi , y j ) = ϕ(xi − Δx, y j ) + Δxg(xi , y j )
(9b)
ϕ(xi , y j ) = ϕ(xi , y j + Δy) + Δyg(xi , y j )
(9c)
ϕ(xi , y j ) = ϕ(xi , y j − Δy) + Δyg(xi , y j )
y
j=N j=N-1
φ1,N-1
φM-1,N-1
j=2
j=1 φ11 φ21
i=0 i=1 i=2
φM-1,1
x
i=M-1 i=M
各节点处的差分方程可由式(7)得到：
[ ] [ ] φi+1, j − 2φij + φi−1, j Δx2 + φi, j+1 − 2φij + φi, j−1 Δy2 = 0
四、实验要求
1、理解并掌握有限差分法的基本原理和处理方法，独立完成程序的编写与调试。 2、完成全部实验内容，并在规定的时间内提交实验报告。
6
A(k,k+1)=1.0/dx^2;
A(k,k+(m-1))=1.0/dy^2;
elseif j==N-1
A(k,k-1)=1.0/dx^2;
A(k,k)=-2.0*(1.0/dx^2+1.0/dy^2);
A(k,k+1)=1.0/dx^2;
A(k,k-(m-1))=1.0/dy^2; b(k)=f(i*dx)/dy^2; end
3、差分方程的求解 如图所示边界条件的情况下确定区域内的静电场电位分布。
y
φ=f(x) b
φ=0ห้องสมุดไป่ตู้
0 0
φ=0
3
φ=0
x a
为了利用有限差分法计算电位的分布，把区域离散成边长为步长 Δx = a M 、 Δy = b N 的长方形网格，如下图。
注意到，在边界 x=0，0<y<b、y=0，0<x< a 以及 x= a，0<y< b 上的电位为常数 0；在边界 y=0，0<x< a 上电位 f（x）也为已知，所以点(i，N) (i=1~M-1)的电位已经 给出，点(i，0)，(i，N)，(0，j)，(M，j)的电位全为零，无需求解。有限差分法的任 务就是确定点(i，j)，(i=1~M-1，j=1~N-1)的电位。
（i=1~M-1，j=1~N-1）
将内部节点按行重新排序：令 k=(j-1)*(M-1)+i， 且 φk = φij ，显然 k=1~(M-1)*(N-1)。
则当 j≠1、N-1，且 i≠1、M-1 时，有
[ ] [ ] φk +1 − 2φk + φk −1 Δx2 + φk +(m−1) − 2φk + φk −(m−1) Δy2 = 0
偏导数 ∂ϕ 可用 ϕ(xi + Δx, y j ) − ϕ(xi , y j ) 来代替，则得
∂x
Δx
ϕ(xi ,
yj
)
− ϕ(xi Δx
+
Δx,
yj
)
=
g(xi ,
yj
)
即
ϕ(xi , y j ) = ϕ(xi + Δx, y j ) + Δxg(xi , y j )
(9a)
类似可得，在边界 x=a、y=0、y=b 上，第二类边界条件的离散形式分别为：
(j=N-1)
写成矩阵形式为
AV = b
（11）
其中系数矩阵 A 为(N-1)*(M-1)阶矩阵，V =[ϕ1,ϕ2,"ϕ(N −1)*(M −1) ]′ 为待求的列向量，
b 是已知的(N-1)*(M-1)列向量。该矩阵方程可由以下循环生成（以 Matlab 语言
为例）。
A=zeros((N-1)*(M-1),(N-1)*(M-1)); % 事先把系数矩阵各元素置0，在重新赋值前保持为0
这就是拉普拉斯方程在节点 (xi , y j ) 处的差分方程。
（7）
2、边界条件的处理
边界条件在电磁场问题的求解中起着至关重要的作用，因为不同电磁场问题可
以满足相同的偏微分方程，它们之所以有不同的解就是因为满足的边界条件不同。
为求得区域内各节点上近似解，还需将原问题中的边界条件化成离散的边界条件，
无源区域静电场的电位满足拉普拉斯方程：∇2ϕ = 0 ，在直角坐标系中，二维表 达式具体为：
∂ 2ϕ ∂x 2
+
∂ 2ϕ ∂y 2
=
0
（1）
有限差分法最关键的一步就是在每个离散点上将微分方程中涉及的各阶导数用
该点的差商近似表示。因为对任意函数 u(x) ，当 Δx 很小时有
u(x + Δx) = u(x) + u′(x)Δx + 1 u′′(x)(Δx)2 + " 2
2
第二类边界条件 ∂ϕ = g(x, y) ：
∂n Γ
以边界
x=0
为例，此时
K n
=
−eKx
，边界上任一节点
(
xi
,
y
j
)
处有
K n
(xi , y j )
∂ϕ ∂n
( xi , y j )
=
⎜⎜⎝⎛
∂ϕ ∂x
K ex
+
∂ϕ ∂y
K ey
⎟⎟⎠⎞
⋅
K n
=
−
∂ϕ ∂x
由方程（2）忽略高次项得： u′(x) ≈ u(x + Δx) − u(x) Δx
(9d)
至于第三类边界条件，在 Γ1和 Γ2 上，可分别仿照式(8)和式(9)处理。
具体求解时，若为第一类边界条件，则边界节点上的近似结果直接由式(8)给出，
不必计算，内部节点近似值由各内部节点的差分方程联立求出；若为第二类边界条
件，则边界节点和内部节点上的近似解由式(7)和(9)联立计算；若为第三类边界条件， 一部分边界节点的近似结果不需计算，直接由式(8)给出，而另一部分边界节点以及 内部节点上的解则需由式(7)和(9)联立求得。这里只要求考虑第一类边值问题。
(10a)
当 i=1 时，有
[ ] [ ] φk +1 − 2φk + 0 Δx2 + φk +(m−1) − 2φk + φk −(m−1) Δy2 = 0
(10b)
当 i= M-1 时，有
[ ] [ ] 0 − 2φk + φk −1 Δx2 + φk +(m−1) − 2φk + φk −(m−1) Δy2 = 0
elseif i==M-1
A(k,k-1)=1.0/dx^2;
A(k,k)=-2.0*(1.0/dx^2+1.0/dy^2);
A(k,k-(m-1))=1.0/dy^2;
A(k,k+(m-1))=1.0/dy^2;
end
elseif j==1
A(k,k-1)=1.0/dx^2;
A(k,k)=-2.0*(1.0/dx^2+1.0/dy^2);