有限差分法及matlab实现

有限差分法解静电场的边值问题的算法实现及相关问题讨论：

王宁远

中国科学技术大学09级物理2班E-mail [email protected]

摘要：

本文用MATLAB实现了有限差分法解静电场边值问题的算法，将偏微分方程的问题化为线性方程组问题，并使用了迭代法进行线性方程组的数值解。讨论了从几个角度去优化迭代法的措施。并运用这样的方法解决了文[1]的闪电模拟问题，，使用了更优化的算法对重新进行了计算，并一定程度上改进了模型,讨论了几个与文[1]所持的不同的观点。

正文：

经典场的边值问题在数学上表达为泊松方程和拉普拉斯方程，但解偏微分方程往往是困难的。幸而很多时候对于具体问题我们需要的不是解析解，而是数值解，所以可以考虑用连续变量离散化的方法求出数值解,在足够的精度上进行逼进,这就引出了有限差分法。

1.1有限差分法:

有限差分法：

微分：f'x=fxh−fxh0=dy

hdx

用有限的h代替，使得

f'x≈

△y△x

差分的种类：

fxh−fxfx−fx−hfxh−fx−h

或者或者

h2hh

fxh−fxfx−fx−h

−

hhfxhfx−h−2fx二阶差分: =

hh2

设Ux,y,z为空间电势的函数。

一阶差分：泊松方程：

∇2U=ρ

使用二阶差分代替泊松方程中的拉普拉斯算符，有：

fxh,y,zfx−h,y,z−2fx,y,z∂2U∂2U∂2U

∇U==∑

∂x2∂y2∂z2h2

2

∑

表示分别对三个变元求差分之和，以下相同

矩阵（数组）是计算机中重要的数据结构，为了方便用矩阵去存储数据，我们网格去划分空间，从而不仅使方程化为简单的有限差分形式，而且这样的数据类型在

计算机中易于储存和运算。那么h=k=l=1，并且令f(x,y,z)=u(x,y,z)就有

Ui1,j,kUi−1,j,kUi,j1,kUi,j−1,kUi,j,k1Ui,j,k−1

=6Ui,j,kρ

这就是泊松方程的有限差分形式，以下估计该方程的精度：由泰勒公式，易知有以下结果：

∂U1∂2U21∂3U3

Uxh,y,z=Ux,y,zhhh⋯23

∂x2∂x6∂x∂U1∂2U21∂3U3

Ux,yk,z=Ux,y,zkkk⋯23

∂y2∂y6∂y

23

∂U1∂U21∂U3

Ux,y,zl=Ux,y,zlll⋯

∂z2∂z26∂z3

若考虑离散的点：

∂U1∂U1∂U

Ui1,j,k=Ui,j,k⋯23

∂x2∂x6∂x∂U1∂2U1∂3U

Ui,j1,k=Ui,j,k⋯

∂y2∂y26∂y3∂U1∂2U1∂3U

Ui,j,k1=Ui,j,k⋯

∂z2∂z26∂z3∂U1∂2U1∂3U

Ui−1,j,k=Ui,j,k−−⋯23

∂x2∂x6∂x∂U1∂2U1∂3U

Ui,j1,k=Ui,j,k−−⋯23

∂y2∂y6∂y∂U1∂2U1∂3U

Ui,j,k1=Ui,j,k−−⋯

∂z2∂z26∂z3

上述六式相加

23

Ui1,j,kUi−1,j,kUi,j1,kUi,j−1,kUi,j,k1Ui,j,k−1

222∂U∂U∂U

⋯=6Ui,j,k222

∂x∂y∂z

2

代入∇U=ρ

Ui1,j,kUi−1,j,kUi,j1,kUi,j−1,kUi,j,k1Ui,j,k−1=6Ui,j,k∇2U⋯=6Ui,j,kρ⋯

由于所有的奇次项被抵消，所以方程：

Ui1,j,kUi−1,j,kUi,j1,kUi,j−1,kUi,j,k1Ui,j,k−1=6Ui,j,kρ ----- (1) 式

的精度为三阶，四阶及更高阶项被略去。 若满足∇U=0有

2

Ui1,j,kUi−1,j,kUi,j1,kUi,j−1,kUi,j,k1Ui,j,k−1

=6Ui,j,k ----- (2)式

此即拉普拉斯方程的有限差分形式。

这里，我们通过有限差分的方法，把偏微分方程在三阶精度下简化为形式易于计算的代数方程。从而使之易于在计算机上实现。

注：有时我们需要解二维静电场，则方程退化为：

∂2U∂2U

Ui1,jUi−1,jUi,j1Ui,j−1=4Ui,j2

∂x∂y2

即

∂2U∂2U

Ui1,jUi−1,jUi,j1Ui,j−1=4Ui,j2

∂x∂y2

1.2算法选择:

下面我们对算法再进行一些讨论。

MATLAB的M语言本身带有矩阵的数据类型，且MATLAB具有高效的数值计算功能。所以我们选择通过MATLAB的M语言去实现这种算法。

(可以看出，三维情况与二维情况没有本质的区别，所以在一下的一系列讨论中，我们为方便起见都只考

虑二维情况。)

由一个简单的例子出发先讨论：

矩形截面由四块导体板围成，其中一块有100v，其他三块全部接地（电势为0），求解平面各处电势：

解：

基于有限差分法，我们得到的差分形式事实上是线性方程组。

由于边界上的点电势已经作为已知条件所固定。所以实际上的未知数只有15×9=135个，而对于每一个未知数（每一个点），我们都可以写出一个方程（每个点的电势等于它旁边四个点的电势的平均值），

从物理意义上考虑，方程确实是独立的，所以方程是可解的，我们那么考虑构造一个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; 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); end

计算出的矩阵：

绘图：

subplot(1,2,1),mesh(v2)axis([0,17,0,11,0,100])

subplot(1,2,2),contour(v2,32)

有限差分法求出的电势的分布图像：

在这个例子中，我们对网格的划分并不细致（17×11），但却要去计算一个135阶方阵的逆，算法的时间复杂度和空间复杂度都比较大，而且对于边界条件较复杂的情况，这样的矩阵会使得程序的编写十分不便，由于我们给出的差分公式本身就是近似公式，求出该方程的精确解没有实际意义，我们只需要一个符合精度要求的近似解即可，我们所以我们应该考虑一个效率更高的，更容易实现的算法，那么可以考虑使用迭代法。

迭代解法:

程序如下：

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 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算法改进

提前使用新值：因为在上述程序中，我们的扫描赋值方式是一行一行扫描，在对v2（i，j）赋值时，它周围四个点其中有2个已经被扫描过了，即已经获得了新的数值，这个数值应该更优，所以我们尽量使用新算出的数值进行迭代，其中：只要把上述代码中

v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v1(i-1,j)+v1(i,j-1))/4改为

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太小时，每次迭代U不能获得足够的增量。而当m太大，则会使得增量过大，在超过目标值时需要更多的迭代次数来返回。那么是否有一种办法能够精确算出最合适的m值或者估计出较合适的m值。

从多次实验看来，当m>=2时计算总是不收敛，而m的最佳取值往往和网格的行列数有关:有资料给出了经验公式:

cosa=

笔者试验,该公式是有效的,所以下文中我们将使用该公式.

同时我们还可以考虑进行扫描方式的优化，上面方法的思路，是尽量使用新算出来的点的电势值.以下我们将使扫描尽量从边界出发，充分利用边界点电势初值的真实性。以下是就这个例子所做的尝试：

当m=1.2迭代次数 逐行扫描：151次 隔行扫描：148 对称扫描：132 当m=1.3迭代次数 逐行扫描122次 隔行扫描：120 对称扫描：104 当m=1.4迭代次数 逐行扫描96次 隔行扫描：94 对称扫描：78次 当m=1.5迭代次数 逐行扫描71次 隔行扫描：69 对称扫描：52次当m=1.6迭代次数 逐行扫描43次 隔行扫描：40 对称扫描：42次当m=1.7迭代次数 逐行扫描62次 隔行扫描：58 对称扫描：75次

可见，隔行扫描能够带来扫描次数的略微减少，而对称扫描在整体看来，一定程度上再次减少了扫描次数，但是在m值较为合适时，并不具有优势。

由上文我们可以知道,这样的算法完全可以在计算机上有效实现,只要给定边界条件,就可以解出拉普拉斯方程或者泊松方程.

我们现在还可以再使用这种算法解决和讨论一个实际的模型：

ππcoslxly;m;21−a

相关问题：

2.1 模型引入:

在《用计算机模拟闪电形成的尝试》--金秀儒 ，（以下称作 文[1] ）一文中，作者试图利用类似的方法解决闪电形成的模型问题。

该文中作者参考了由Niemeyer Pietronero和Wielsmann于1983年提出的介电击穿模型

(dielectric breakdown model，DBM)，结合该问题该模型可简述如下：

该模型全称 Dielectric Breakdown Moel ，是一个从概率角度去描述介质电击穿的模型。正如上图,每一个点都有可能击穿离它最近的四个点,每击穿一个点,击穿点的电势就可以认为和击穿源相同.由于雷雨天气中云层顶部带正电，底部带负电，所以会使得地面成为高电势，云层成为低电势，所以在计算中我们不妨假设云层的电势为0,地面的电势为100.天地间简化为一个平面。闪电的路径为击穿的点，由于闪电从云端一直到地面，第一个被击穿的点在云端，可以认为所有击穿的点的电势都为0,，其概率满足

Pi=

φi

n1

η

-------- (3)式

∑φηi

我们假设第一个击穿的点在云端的中间，每击穿一个点，我们把已被击穿的点计入边界，这满足拉普拉斯方程，所以我们再用有限差分法迭代重新计算平面的电势分布，然后依照新算出来的电势再配合随机函数得到下一个击穿的点。每个已被击穿的点都可能会作为生长点，使最近的四个点之一被击穿。最终击穿点到达边界则计算结束。

因为计算电势需要完整的边界条件，所以我们做一个近似的假定，我们认为在整个过程中平面左右两个边界上的电势一直不会改变，即满足最初的线性变化的分布。

2.2 模型讨论与改进:

在这里我们利用迭代的算法完成作者未能完成的设计和计算，并且指出文中的几点缺陷。

首先文[1]认为

“最自然的想法是直接通过递归的方法用Average函数求解，但实现过程中发现，这个方法是完全行不通的，因为它会产生所谓循环引用的问题，

即若干步后所求变量计算时调用自己的值，使程序进入死循环。”

“以及由于方程中多为0元，为稀疏方程组，故可以用雅可比迭代法计算，但由于要判断迭代是否收敛以及收敛速度的问题，实现时可能需要用并行算法，皆超出作

者能力范围，而且计算机的性能上也不符合条件，需要进一步优化解法，需要加强计算机功能；”

这是不妥的，因为在这里实在没有必要运用递归求解，而我们正是利用迭代法解决了这个问题，而且前文中已经验证，在解决拉普拉斯方程的有限差分形式的问题中，迭代算法比直接解方程组的算法拥有更低的时间复杂度。文[1]认为，在

“η=1时，模型满足Laplace方程”

我们认为这也是不妥的，因为首先作者没有写出这样做的理论依据，而且在最初提出该模型的文[5]中，也没有指出必须只考虑η =1的情况，而且我们将会在下文指出，当η>1时，尤其是 η在2

左右

时，我们能得到

更为接近真实情况的解。

在文[1]中，作者最大只实现了20×20像素的情况，离模拟真实情况相距甚远，这里我们运用更优的算法（超松弛迭代法）能够计算更高像素的解，将会更接近真实情况。

同时，再考察DBM模型中的核心公式(3)式，由于对介质电击穿的具体原理以及性质不甚了解，所以DBM模型选择通过电势计算概率来描述这样的现象，但是笔者认为，更确切的描述应该是：

Pi=

Ei

n1

η

-

∑Eηi

进一步说,击穿本来就是对稠密的点的连续击穿过程,我们写成Ei应当只有在对于足够近的两个点计算时才是精确的，不过我们这里假定网格的划分已经足够细致，因而这样的近似的是完全合理的。同时因为我们方程中处理的都是电势,所以改写(4)试可以得到:

Pi=

φi/di

n1

η

-------- (4)式

∑φi/diη

很明显，当网格均匀时，所有di都相等，(4)式就蜕变为(3)式。所以，根据(4)式我们可以重新考察击穿的形式

在原先的模型中，从点(i,j)出发，仅仅考虑了对(i+1,j),(i-1,j),(i,j-1),(i,j+1)点的击穿，我们用(4)式代替(3)式，再考虑加入(i+1,j-1),(i-1,j-1),(i+1,j+1),(i-,j+1)四个点，考虑到(4)式中距离di的影响，我们对这四个点只要把φ改写为φ/就可以直接应用原来的式(3)求解。相对于文[1]作者所描述，对于20×20以上的情况：

“由于电脑功能尚不够强大，运行容易导致系统崩溃，无法实现。”于是：

在笔者的计算机软硬件环境下,进行了测试： CPU：Core2 Duo 2.1Ghz OS：Ubuntu10.04 64bit Soft：Matlab2009a for Unix用时计算：

η=3

可以看出，我们不仅轻松实现了20×20以上的情况，而且最高实现了300×300的情况,将像素从400提高到90000,我们的运算效率有着质的改变。

这里笔者猜想，文[1]作者由于采取了直接解方程组的办法，那么对于100×100的空间，编需要一个10000×10000的浮点系数矩阵，这将占据非常大的内存空间，所以文[1]所说的容易导致系统崩溃也就不难解释了。而且在文[1]中，作者提出了如下遗憾： ·程序未涉及电势变化

借助于有限差分法，这里我们解决了这个问题。

2.3 绘图与进一步讨论:

这里二维和三维没有本质的区别，简洁起见我们只考虑二维的情形，我们使用MATLAB的M语言编写程序，代码放在最后。

下面我们绘出我们所算出的图形（该图形先由MATLAB依矩阵元素取值绘出，然后用GIMP上色得到），并讨论η的取值对最终结果的影响：

插图

1: 300*300;η=3

插图 2: : 250*250;η=2

插图 3: 200*200;η=1

我们这里可以看出来， 随着η

的增大，闪电的“分叉”情况越来越少，根据生活经验， η大约取2或3之间时应该更符合客观情况。

基

于这样的结果，我们还可以讨论一下地面上的物体对闪电轨迹的影响，例如，加入以下代码：

v(70:100,70)=1;

v1(70:100,70)=100;

trace(70:100,70)=2;

重写标记矩阵，在地上（70:100,70）插上一根旗杆（避雷针）。

我们多次运行程序，闪电基本上都能能打到旗杆上，说明地面高耸尖端的物体的确会吸引闪电,或者说”上帝喜欢把任何灾难都降落在那些自高自大的东西上面.”

由于时间和精力原因，这里没有做更多的讨论，事实上我们还可以讨论不同形状的物体的引雷区域，引雷概率等，以及该有限差分方法在更多模型上的应用等等。

2.4 附录：

附:程序,由M语言编写.

tic;

lx=100;ly=100;

a=(cos(pi/lx)+cos(pi/ly))/2;

w=2/(1+sqrt(1-a*a));

v1=zeros(ly,lx);

v=v1;

r=3;

trace=v;

for i=1:ly

for j=1:lx

v1(i,j)=(i-1)*100/(ly-1);

end

end

o=sqrt(2);

v(1,:)=1;v(ly,:)=1;

v(:,1)=1;v(:,lx)=1;

k=0;

i=2;j=fix(lx/2);

while(v(i,j)~=1)

v(i,j)=2; trace(i,j)=1;

v1(i,j)=0;

v2=v1;maxt=1;

while(maxt>1e-5)

k=k+1;

maxt=0;

for m=2:ly-1

for n=2:lx-1;

if v(m,n)==0

v2(m,n)=v2(m,n)+w*(v1(m,n+1)+v1(m+1,n)+v2(m-1,n)+v2(m,n-1)-4*v2(m,n))/4; t=abs(v2(m,n)-v1(m,n));

if(t>maxt) maxt=t;

end

end

end

end

v1=v2;

end

p=0;

for i=2:ly-1

for j=2:lx-1

if (v(i,j)==2)

if ( v(i-1,j)==2) p1=0; else p1=v1(i-1,j)^r; end;

if (v(i+1,j)==2) p2=0; else p2=v1(i+1,j)^r;end;

if (v(i,j-1)==2) p3=0; else p3=v1(i,j-1)^r;end;

if (v(i,j+1)==2) p4=0; else p4=v1(i,j+1)^r;end;

if ( v(i-1,j-1)==2) p5=0; else p5=(v1(i-1,j-1)/o)^r; end;

if ( v(i-1,j+1)==2) p6=0; else p6=(v1(i-1,j+1)/o)^r; end;

if ( v(i+1,j-1)==2) p7=0; else p7=(v1(i+1,j-1)/o)^r; end;

if ( v(i+1,j+1)==2) p8=0; else p8=(v1(i+1,j+1)/o)^r; end;

p=p1+p2+p3+p4+p5+p6+p7+p8+p;

end

end;

end;

q=rand(1,1);

s=0;

for i=2:ly-1

for j=2:lx-1

if (v(i,j)==2)

if ( v(i-1,j)==2) p1=0; else p1=(v1(i-1,j)^r/p); end;

if (v(i+1,j)==2) p2=0; else p2=(v1(i+1,j)^r/p);end;

if (v(i,j-1)==2) p3=0; else p3=(v1(i,j-1)^r/p);end;

if (v(i,j+1)==2) p4=0; else p4=(v1(i,j+1)^r/p);end;

if ( v(i-1,j-1)==2) p5=0; else p5=((v1(i-1,j-1)/o)^r)/p; end;

if ( v(i-1,j+1)==2) p6=0; else p6=((v1(i-1,j+1)/o)^r)/p; end;

if ( v(i+1,j-1)==2) p7=0; else p7=((v1(i+1,j-1)/o)^r)/p; end;

if ( v(i+1,j+1)==2) p8=0; else p8=((v1(i+1,j+1)/o)^r)/p; end;

s=s+p1+p2+p3+p4+p5+p6+p7+p8;

if (q>=s-p1-p2-p3-p4-p5-p6-p7-p8) && (q

m=i-1; n=j;

end;

if (q>=s-p2-p3-p4-p5-p6-p7-p8) && (q

m=i+1;n=j;

end;

if (q>=s-p3-p4-p5-p6-p7-p8) && (q

m=i;n=j-1;

end;

if (q>=s-p4-p5-p6-p7-p8) && (q

m=i;n=j+1;

end;

if (q>=s-p5-p6-p7-p8) && (q

m=i-1;n=j-1;

end;

if (q>=s-p6-p7-p8) && (q

m=i-1;n=j+1;

end;

if (q>=s-p7-p8) && (q

m=i+1;n=j-1;

end;

if (q>=s-p8) && (q

m=i+1;n=j+1;

end;

end

end

end

i=m;j=n;

end

k

toc;

time=toc;

im=pcolor(trace)

colormap(cool);set(im, 'LineStyle','none');

参考文献：

[1] 《用计算机模拟闪电形成的尝试》 金秀儒 2006

[2] 《数值计算方法与算法》 科学出版社 张韵华 奚梅成 陈效群 编

[3] 《偏微分数值解法基础教程》科学出版社 林群 编著

[4] 《数值解法和MATLAB实现与应用》Gerald Recktenwald 著 伍卫国 万群 张辉 等译

[5] Fractal Dimension of Dielectric Breakdown Niemeyer Pietronero

Wielsmann 1983

[6] MATLAB 6 实例教程 郝红伟 编著 中国电力出版社

[7] 3D modelling of Electrical Discharges D. I. Amarasinghe and D. U. J. Sonnadara 2007

后记：

这篇论文的写作从5月20日开始，一直写到交论文期限前一天的深夜，开始是在数值计算类的书上看到了微分方程离散化的方法，然后试着将这种方法在计算机上实现，后来一直在思考如何应用这种算法做出一些有意义的内容。开始计划了好几个题目，但后来都被自己否定。

直到后来看到了学长的论文，学长模拟的闪电虽然简陋，却给我以十分深刻的印象。一是因为想不到闪电这样的“复杂”现象的模型居然可以这么简洁;二是因为这个模型正好可以利用我的算法进行处理。我觉得自己可以比他做的更好，于是便着手进行了编程。之所以选择MATLAB，是因为它由着较强大的数值计算功能。

做到最后，再回过头来看这篇论文，虽然事实上少有创新的东西。但这三个星期，给了我很大的启发，为了能够理解各种数学方法，到图书馆借了数本相关的书，颇费一番功夫才弄懂的，为了理解DBM模型，还需查阅国外的论文（国内这方面的论文几乎找不到）。同时MATLAB也是为了写这篇论文才进行了较深入的学习。很多东西是需要的应用中学习的，而不单单拘泥于课本。