有限差分法及matlab实现

⋯
U
i,
j
, k1=U
i
,
j
, k
∂U ∂z
12
∂2 U ∂ z2
1 6
∂3U ∂ z3
⋯
U
i−1,
j , k=U i ,
j
,
k−
∂U ∂x
1 2
∂2 U ∂ x2
− 16
∂3 U ∂ x3
⋯
U
i,
j1,k=U i ,
j
,
k−
∂U ∂y
1 2
∂2 U ∂ y2
− 16
∂3 U ∂ y3
⋯
U i ,
f
xh− h
f
x
− h
f
x
−
f h
x−h
=
f
xh
f
x−h−2f h2
x
设 U x , y , z为空间电势的函数 。
泊松方程：
∇2U=ρ
使用二阶差分代替泊松方程中的拉普拉斯算符，有：
∑ ∇
2
U
=
∂2 U ∂ x2
∂2 U ∂ y2
Baidu Nhomakorabea
∂2 U ∂ z2
=
f xh, y , z f x−h , y , z−2f x , y , z h2
改为
v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v2(i-1,j)+v2(i,j-1))/4 解释执行后，k=222;
迭代次数接近原来的一半，可见这样的算法优化是有效的。
为了能使迭代次数进一步减少，考虑让每次迭代获得更多的增量，
那么将：v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v2(i-1,j)+v2(i,j-1))/4; 变形为：v2(i,j)=v1(i,j)+(v1(i,j+1)+v1(i+1,j)+v2(i-1,j)+v2(i,j-1)-4*v1(i,j))/4; 再考虑：v2(i,j)=v1(i,j)+m*(v1(i,j+1)+v1(i+1,j)+v2(i-1,j)+v2(i,j-1)-4*v1(i,j))/4 适当改变 m 的值是否能够减少迭代次数？ 我们做了如下试验:
而当 m 太大，则会使得增量过大，在超过目标值时需要更多的迭代次数来返回。
那么是否有一种办法能够精确算出最合适的 m 值或者估计出较合适的 m 值。
从多次实验看来，当 m>=2 时计算总是不收敛，而 m 的最佳取值往往和网格的行列数有关:
有资料给出了经验公式:
a=
cos
π lx
cos 2
π ly
;
样的方法解决了文[1]的闪电模拟问题，，使用了更优化的算法对重新进行了计算，并一定程度上改进了 模型,讨论了几个与文[1]所持的不同的观点。
正文：
经典场的边值问题在数学上表达为泊松方程和拉普拉斯方程，但解偏微分方程往往是困难的。幸而
很多时候对于具体问题我们需要的不是解析解，而是数值解，所以可以考虑用连续变量离散化的方法求
从物理意义上考虑，方程确实是独立的，所以方程是可解的，我们那么考虑构造一个 135 阶的方阵，只
要求出其逆就可以了。
传统解法:
程序如下：
a=zeros(135,135); for i=1:135 a(i,i)=1;end; for i=1:7 a(15*i+1,15*i+2)=-0.25;a(15*i+1,15*i+16)=-0.25;a(15*i+1,15*i-14)=-0.25; end for i=1:7 a(15*i+15,15*i+14)=-0.25;a(15*i+15,15*i+30)=-0.25;a(15*i+15,15*i)=-0.25; end a(1,2)=-0.25;a(1,16)=-0.25; a(121,122)=-0.25;a(121,106)=-0.25; a(135,134)=-0.25;a(135,120)=-0.25; a(15,14)=-0.25;a(15,30)=-0.25;
for i=2:14 a(i,i-1)=-0.25;a(i,i+1)=-0.25; a(i,i+15)=-0.25; end for i=122:134 a(i,i-1)=-0.25;a(i,i+1)=-0.25; a(i,i-15)=-0.25; end for i=1:7
for j=2:14; a(15*i+j,15*i+j-1)=-0.25; a(15*i+j,15*i+j+1)=-0.25; a(15*i+j,15*i+j+15)=-0.25; a(15*i+j,15*i+j-15)=-0.25;
for j=2:lx-1;
v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v1(i-1,j)+v1(i,j-1))/4; %进行迭代计算
t=abs(v2(i,j)-v1(i,j));
if(t>maxt) maxt=t;end
end
end
v1=v2;
end k=419
%输出迭代次数
1.3 算法改进
这里，我们通过有限差分的方法，把偏微分方程在三阶精度下简化为形式易于计算的代数方程。从而使 之易于在计算机上实现。
注：有时我们需要解二维静电场，则方程退化为：
U i1,
jU i−1,
jU i ,
j1U i ,
j−1=4Ui ,
j
∂2U ∂ x2
∂2 U ∂ y2
即
U
i1,
jU i−1,
jU i ,
lx=17;ly=11;
%定义矩阵维数
v1=zeros(ly,lx); for j=2:lx-1 v1(ly,j)=100;
%建立一个矩阵
end
v2=v1;maxt=1;t=0; k=0;
%设置边界条件
while(maxt>1e-6)
%精度要求，达到精度要求跳出循环
k=k+1
maxt=0;
for i=2:ly-1
k3⋯
U
x
,
y
,
z l =U
x
,
y
,
z
∂U ∂z
l
1 2
∂2 U ∂ z2
l2
1 6
∂3U ∂ z3
l3⋯
若考虑离散的点：
U
i1,
j , k=U i
,
j
, k
∂U ∂x
1 2
∂2 U ∂ x2
1 6
∂3U ∂ x3
⋯
U
i,
j1, k=U
i
,
j
, k
∂U ∂y
1 2
∂2 U ∂ y2
1 6
∂3U ∂ y3
∑ 表示分别对三个变元求差分之和，以下相同
矩阵（数组）是计算机中重要的数据结构，为了方便用矩阵去存储数据，我们网格去划分空间，从而不 仅使方程化为简单的有限差分形式，而且这样的数据类型在
计算机中易于储存和运算。那么 h=k=l=1，并且令 f(x,y,z)=u(x,y,z)
就有
U i1, j , kU i−1, j , kU i , j1, kU i , j−1,kU i , j , k1U i , j , k−1
出数值解,在足够的精度上进行逼进,这就引出了有限差分法。
1.1 有限差分法:
有限差分法：
微分：
f
' x=
f
xh− h
f
x
h 0=
dy dx
用有限的 h 代替，使得
f
'
x≈
△ △
y x
差分的种类：
一阶差分：
f xh− f x h
或者
f x −f x−h h
或者
f xh− f x−h 2h
二阶差分:
∂2U ∂ z2
⋯
代入 ∇2 U =ρ
U i1, j , kU i−1, j , kU i , j1, kU i , j−1,kU i , j , k1U i , j , k−1
= 6Ui , j , k∇2 U ⋯=6Ui , j ,kρ⋯
由于所有的奇次项被抵消，所以方程：
U i1, j , kU i−1, j , kUi , j1, kU i , j−1, kU i , j , k1U i , j , k−1
m
=
2
1 1 −a2
;
笔者试验,该公式是有效的,所以下文中我们将使用该公式.
同时我们还可以考虑进行扫描方式的优化，上面方法的思路，是尽量使用新算出来的点的电势值.以下
我们将使扫描尽量从边界出发，充分利用边界点电势初值的真实性。
以下是就这个例子所做的尝试：
当 m=1.2 迭代次数 逐行扫描：151 次 隔行扫描：148 对称扫描：132
有限差分法解静电场的边值问题的算法实现及相关问题讨论：
王宁远
中国科学技术大学 09 级物理 2 班 E-mail wny@mail.ustc.edu.cn
摘要：
本文用 MATLAB 实现了有限差分法解静电场边值问题的算法，将偏微分方程的问题化为线性方程
组问题，并使用了迭代法进行线性方程组的数值解。讨论了从几个角度去优化迭代法的措施。并运用这
j , k1=U i ,
j
,
k−
∂U ∂z
12
∂2 U ∂ z2
−
1 6
∂3 U ∂ z3
⋯
上述六式相加
U i1, j , kU i−1, j , kUi , j1, kU i , j−1, kU i , j , k1U i , j , k−1
=
6Ui ,
j
,
k
∂2 U ∂ x2
∂2 U ∂ y2
m
1.2
k(迭代
次数）
151 次
1.3 122 次
1.4 96 次
1.5 71 次
1.6 43 次
1.7 62 次
1.8 88 次
>2
不收敛
可见，这样的更改在 m 取合适的值的时候能带来迭代次数十分显著的减少，在数值计算中,这样的算法称
作超松弛迭代算法.
但什么样的 m 才是“合适的”值，因为当 m 太小时，每次迭代 U 不能获得足够的增量。
j1U i ,
j−1=4Ui ,
j
∂2U ∂ x2
∂2 U ∂ y2
1.2 算法选择:
下面我们对算法再进行一些讨论。
MATLAB 的 M 语言本身带有矩阵的数据类型，且 MATLAB 具有高效的数值计算功能。所以我 们选择通过 MATLAB 的 M 语言去实现这种算法。
(可以看出，三维情况与二维情况没有本质的区别，所以在一下的一系列讨论中，我们为方便起见都只考
时间复杂度和空间复杂度都比较大，而且对于边界条件较复杂的情况，这样的矩阵会使得程序的编写十
分不便，由于我们给出的差分公式本身就是近似公式，求出该方程的精确解没有实际意义，我们只需要
一个符合精度要求的近似解即可，我们所以我们应该考虑一个效率更高的，更容易实现的算法，那么可
以考虑使用迭代法。
迭代解法:
程序如下：
虑二维情况。)
由一个简单的例子出发先讨论：
矩形截面由四块导体板围成，其中一块有 100v，其他三块全部接地（电势为 0），求解平面
各处电势：
解： 基于有限差分法，我们得到的差分形式事实上是线性方程组。
由于边界上的点电势已经作为已知条件所固定。所以实际上的未知数只有 15×9=135 个，而对于
每一个未知数（每一个点），我们都可以写出一个方程（每个点的电势等于它旁边四个点的电势的平均 值），
提前使用新值：因为在上述程序中，我们的扫描赋值方式是一行一行扫描，在对 v2（i，j）赋值时，它 周围四个点其中有 2 个已经被扫描过了，即已经获得了新的数值，这个数值应该更优，所以我们尽量使
用新算出的数值进行迭代，其中：
只要把上述代码中
v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v1(i-1,j)+v1(i,j-1))/4
= 6Ui , j , kρ
这就是泊松方程的有限差分形式，以下估计该方程的精度： 由泰勒公式，易知有以下结果：
U
xh
,
y
,
z =U
x
,
y
,
z
∂U ∂x
h
1 2
∂2U ∂ x2
h2
1 6
∂3 U ∂ x3
h3⋯
U
x
,
yk,
z =U
x
,
y
,
z
∂U ∂y
k
1 2
∂2U ∂ y2
k2
1 6
∂3 U ∂ y3
= 6Ui , j , kρ
----- (1) 式
的精度为三阶，四阶及更高阶项被略去。
若满足 ∇2 U =0
有
U i1, j , kU i−1, j , kU i , j1, kU i , j−1,kU i , j , k1U i , j , k−1
= 6Ui , j , k
----- (2)式
此即拉普拉斯方程的有限差分形式。
当 m=1.3 迭代次数 逐行扫描 122 次 隔行扫描：120 对称扫描：104
当 m=1.4 迭代次数 逐行扫描 96 次 隔行扫描：94 对称扫描：78 次
当 m=1.5 迭代次数 逐行扫描 71 次 隔行扫描：69 对称扫描：52 次
当 m=1.6 迭代次数 逐行扫描 43 次 隔行扫描：40 对称扫描：42 次
end
计算出的矩阵：
绘图：
subplot(1,2,1),mesh(v2) axis([0,17,0,11,0,100]) subplot(1,2,2),contour(v2,32)
有限差分法求出的电势的分布图像：
在这个例子中，我们对网格的划分并不细致（17×11），但却要去计算一个 135 阶方阵的逆，算法的
end end b=a^(-1); c=zeros(135,1); for i=121:135 c(i,1)=25;end d=b*c;s=zeros(11,17);for i=2:16 s(11,i)=100; end for i=1:9
for j=1:15; s(i+1,j+1)=d(15*(i-1)+j,1);