基于Matlab和高斯-赛德尔迭代应用有限差分法求电势及电场分布
Posted 迟钝皮纳德
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了基于Matlab和高斯-赛德尔迭代应用有限差分法求电势及电场分布相关的知识，希望对你有一定的参考价值。
基于Matlab和高斯-赛德尔迭代应用有限差分法求电势及电场分布
电磁场的边界条件
根据电磁场理论，若想知道某一区域内的电场及电势分布，我们只需要知道其边界条件处的电势情况就可以了。这样的边界条件可以使直接给出边界处的电势值，也可以给出电势在某一方向的方向导数的值或是其他的面积分的形式。但是这其中比较关键的一点就是需要求得某一封闭区域内的电场，那么这封闭区域边界的所有电势情况必须全部了解，而不能有疏漏。
我们考虑二维的情况，对于静电场的边值条件：
拉普拉斯方程
▽
2
φ
=
∂
2
φ
∂
x
2
+
∂
2
φ
∂
y
2
=
−
ρ
ε
\\bigtriangledown^2\\varphi = \\frac\\partial^2\\varphi\\partial x^2 + \\frac\\partial^2\\varphi\\partial y^2 = -\\frac\\rho \\varepsilon
▽2φ=∂x2∂2φ+∂y2∂2φ=−ερ
φ
∣
S
=
f
1
(
s
)
\\varphi |_S =f_1 (s)
φ∣S=f1(s)
∂
φ
∂
n
⃗
∣
S
=
f
2
(
s
)
\\frac\\partial \\varphi \\partial \\vecn |_S=f_2 (s)
∂n∂φ∣S=f2(s)
或者是线性组合
(
α
φ
+
β
∂
φ
∂
n
⃗
)
∣
S
=
f
3
(
s
)
(\\alpha \\varphi +\\beta \\frac\\partial \\varphi \\partial \\vecn )|_S =f_3 (s)
(αφ+β∂n∂φ)∣S=f3(s)
如果上述条件在边界处给完全确定，那么边界内部的任何一点电势也就完全确定了。但是这并不排除某些不可解的情况，例如边界区域电势为常数。这种情况下，通过该算法所确定的解有无穷多个，因为这样边界条件下，内部电势既可以为常数，或是包含某些奇点的情况，例如中心有一个点电荷。
以上情况我们假定的是该边界内部介质条件不变。如果条件不是单一的，例如中心存在点电荷或电介质变化。我们需要对该区域分割成多个区域来分别求解内部电势分布。
高斯-赛德尔迭代算法
我们考虑二维的情况。由于电势为标量。为了保证场内区域电势的连续性，对于高斯-赛德尔迭代算法的理解可以近似成该点的电势值，可以由周围几个点电势的平均值来代替。
φ
i
,
j
(
k
+
1
)
=
1
4
(
φ
i
−
1
,
j
(
k
)
+
φ
i
+
1
,
j
(
k
)
+
φ
i
,
j
−
1
(
k
)
+
φ
i
,
j
+
1
(
k
)
)
\\varphi _i,j^(k+1) =\\frac14 (\\varphi _i-1,j^(k) +\\varphi _i+1,j^(k)+\\varphi _i,j-1^(k)+\\varphi _i,j+1^(k))
φi,j(k+1)=41(φi−1,j(k)+φi+1,j(k)+φi,j−1(k)+φi,j+1(k))
k
k
k表示迭代次数。
将场域内进行方格坐标化之后，对场域内每个点进行迭代。同时保持边界处的值不变，最后场内任意点在某次迭代之后数值会趋于收敛。即满足
∣
φ
i
,
j
(
k
+
1
)
−
φ
i
,
j
(
k
)
∣
<
ε
|\\varphi _i,j^(k+1) -\\varphi _i,j^(k) |<\\varepsilon
∣φi,j(k+1)−φi,j(k)∣<ε
对任意
i
,
j
i,j
i,j及
ε
>
0
\\varepsilon>0
ε>0总存在
k
k
k使得该式成立，这里
k
k
k就是迭代次数。
算法举例
内部存在一个矩形区域等势体
本项目为一矩形区域，中心矩形区域导体电势为 U 0 U_0 U0，这里我们设为 100 100 100.导体上为等势体，外围矩形区域电势为 0 0 0。通过高斯-赛德尔迭代法求得中间电势分布，下面可以直接改参数值。由于需要方便矩阵运算，坐标值需要整数。因此若要形成某一区域的电位需要先进行线性变换在求解，最后还可以进行线性插值。
在求得电势后，可以对其求梯度便形成了电场线
E
⃗
=
−
▽
φ
\\vecE =-\\bigtriangledown \\varphi
E=−▽φ
下面的代码可以生成电势图和电场线：
%%
clc;clear;
U0=100; %设置中心区域电位
err=0.001; %设置误差
boolsta=0; %设置判断变量，判断是否停止迭代
[X1,X2,Y1,Y2]=deal(1,30,1,30); %设置外边界区域，似乎只能设置正方形
[x1,x2,y1,y2]=deal(10,20,20,25); %设置内边界区域
itera=0; %迭代次数
A=zeros(X2-X1+1,Y2-Y1+1); %初始化矩阵
A(x1:x2,y1:y2)=U0*ones(x2-x1+1,y2-y1+1); %设置中心区域电位
A1=A; %每次迭代更新后的矩阵
while ~boolsta
for i=X1+1:X2-1
for j=Y1+1:Y2-1
%高斯-塞德尔迭代
if i>=x1 & i<=x2 & j>=y1 & j<=y2
A1(i,j) = U0;
else
A1(i,j) = ( A(i+1,j) + A(i-1,j) + A(i,j+1) + A(i,j-1))/4;
end
%判断是否所有点都在误差范围内
if abs(A1(i,j)-A(i,j))<err
boolsta = boolsta+1;
end
end
end
boolsta = floor( boolsta/(X2-X1-1)/(Y2-Y1-1));
A=A1;
itera = itera +1;
end
itera %显示迭代次数
subplot(1,2,1);
[y,x]=meshgrid(X1:X2 , Y1:Y2);
surf(x,y,A);
xlabel('x');
ylabel('y');
zlabel('电势U');
title('迭代后的电势图');
subplot(1,2,2);
plot([x1,x2,x2,x1,x1],[y1,y1,y2,y2,y1],'k',...
[X1,X2,X2,X1,X1],[Y1,Y1,Y2,Y2,Y1],'k');
xlabel('x');
ylabel('y');
hold on;
h=contour(x,y,A,20);
clabel(h);
text((x1+x2)/2,(y1+y2)/2,'内边界');
text(0.9*X2-0.1*X1,(Y1+Y2)/2,'外边界');
title('等位线图&电场线图');
[dy,dx]=gradient(A);
hold on;
quiver(x,y,-dx,-dy);
输出结果为：
迭代次数
i
t
e
r
a
=
661
itera =661
itera=661次。
内部存在两个（多个）矩形区域等势体
对上面的代码稍微改动，修改初始条件如下：
[U01,U02]=deal(-50,100); %设置中心区域1,2电位
err=0.001; %设置误差
boolsta=0; %设置判断变量，判断是否停止迭代
[X1,X2,Y1,《数值分析》-- 雅可比迭代法高斯—塞德尔迭代法