有限差分法及matlab实现

有限差分法解静电场的边值问题的算法实现及相关问题讨论：王宁远中国科学技术大学09级物理2班[e-mail]摘要：本文用MATLAB实现了有限差分法解静电场边值问题的算法，将偏微分方程的问题化为线性方程组问题，并使用了迭代法进行线性方程组的数值解。讨论了从几个角度去优化迭代法的措施。并运用这样的方法解决了文[1]的闪电模拟问题，，使用了更优化的算法对重新进行了计算，并一定程度上改进了模型,讨论了几个与文[1]所持的不同的观点。正文：经典场的边值问题在数学上表达为泊松方程和拉普拉斯方程，但解偏微分方程往往是困难的。幸而很多时候对于具体问题我们需要的不是解析解，而是数值解，所以可以考虑用连续变量离散化的方法求出数值解,在足够的精度上进行逼进,这就引出了有限差分法。1.1有限差分法:有限差分法：微分：f'x=fxh−fxhh0=dydx用有限的h代替，使得f'x≈△y△x差分的种类：一阶差分：fxh−fxh或者fx−fx−hh或者fxh−fx−h2h二阶差分:fxh−fxh−fx−fx−hhh=fxhfx−h−2fxh2设Ux,y,z为空间电势的函数。泊松方程：∇2U=ρ使用二阶差分代替泊松方程中的拉普拉斯算符，有：∇2U=∂2U∂x2∂2U∂y2∂2U∂z2=∑fxh,y,zfx−h,y,z−2fx,y,zh2∑表示分别对三个变元求差分之和，以下相同矩阵（数组）是计算机中重要的数据结构，为了方便用矩阵去存储数据，我们网格去划分空间，从而不仅使方程化为简单的有限差分形式，而且这样的数据类型在计算机中易于储存和运算。那么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ρ这就是泊松方程的有限差分形式，以下估计该方程的精度：由泰勒公式，易知有以下结果：Uxh,y,z=Ux,y,z∂U∂xh12∂2U∂x2h216∂3U∂x3h3⋯Ux,yk,z=Ux,y,z∂U∂yk12∂2U∂y2k216∂3U∂y3k3⋯Ux,y,zl=Ux,y,z∂U∂zl12∂2U∂z2l216∂3U∂z3l3⋯若考虑离散的点：Ui1,j,k=Ui,j,k∂U∂x12∂2U∂x216∂3U∂x3⋯Ui,j1,k=Ui,j,k∂U∂y12∂2U∂y216∂3U∂y3⋯Ui,j,k1=Ui,j,k∂U∂z12∂2U∂z216∂3U∂z3⋯Ui−1,j,k=Ui,j,k−∂U∂x12∂2U∂x2−16∂3U∂x3⋯Ui,j1,k=Ui,j,k−∂U∂y12∂2U∂y2−16∂3U∂y3⋯Ui,j,k1=Ui,j,k−∂U∂z12∂2U∂z2−16∂3U∂z3⋯上述六式相加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∂x2∂2U∂y2∂2U∂z2⋯代入∇2U=ρ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)式的精度为三阶，四阶及更高阶项被略去。若满足∇2U=0有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)式此即拉普拉斯方程的有限差分形式。这里，我们通过有限差分的方法，把偏微分方程在三阶精度下简化为形式易于计算的代数方程。从而使之易于在计算机上实现。注：有时我们需要解二维静电场，则方程退化为：Ui1,jUi−1,jUi,j1Ui,j−1=4Ui,j∂2U∂x2∂2U∂y2即Ui1,jUi−1,jUi,j1Ui,j−1=4Ui,j∂2U∂x2∂2U∂y21.2算法选择:下面我们对算法再进行一些讨论。MATLAB的M语言本身带有矩阵的数据类型，且MATLAB具有高效的数值计算功能。所以我们选择通过MATLAB的M语言去实现这种算法。(可以看出，三维情况与二维情况没有本质的区别，所以在一下的一系列讨论中，我们为方便起见都只考虑二维情况。)由一个简单的例子出发先讨论：矩形截面由四块导体板围成，其中一块有100v，其他三块全部接地（电势为0），求解平面各处电势：解：基于有限差分法，我们得到的差分形式事实上是线性方程组。由于边界上的点电势已经作为已知条件所固定。所以实际上的未知数只有15×9=135个，而对于每一个未知数（每一个点），我们都可以写出一个方程（每个点的电势等于它旁边四个点的电势的平均值），从物理意义上考虑，方程确实是独立的，所以方程是可解的，我们那么考虑构造一个135阶的方阵，只要求出其逆就可以了。传统解法:程序如下：a=zeros(135,135);fori=1:135a(i,i)=1;end;fori=1:7a(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;endfori=1:7a(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;enda(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;fori=2:14a(i,i-1)=-0.25;a(i,i+1)=-0.25;a(i,i+15)=-0.25;endfori=122:134a(i,i-1)=-0.25;a(i,i+1)=-0.25;a(i,i-15)=-0.25;endfori=1:7forj=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;endendb=a^(-1);c=zeros(135,1);fori=121:135c(i,1)=25;endd=b*c;s=zeros(11,17);fori=2:16s(11,i)=100;endfori=1:9forj=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);%建立一个矩阵forj=2:lx-1v1(ly,j)=100;end%设置边界条件v2=v1;maxt=1;t=0;k=0;while(maxt1e-6)%精度要求，达到精度要求跳出循环k=k+1maxt=0;fori=2:ly-1forj=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(tmaxt)maxt=t;endendendv1=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的值是否能够减少迭代次数？我们做了如下试验:m1.21.31.41.51.61.71.82k(迭代次数）151次122次96次71次43次62次88次不收敛可见，这样的更改在m取合适的值的时候能带来迭代次数十分显著的减少，在数值计算中,这样的算法称作超松弛迭代算法.但什么样的m才是“合适的”值，因为当m太小时，每次迭代U不能获得足够的增量。而当m太大，则会使得增量过大，在超过目标值时需要更多的迭代次数来返回。那么是否有一种办法能够精确算出最合适的m值或者估计出较合适的m值。从多次实验看来，当m=2时计算总是不收敛，而m的最佳取值往往和网格的行列数有关:有资料给出了经验公式:a=cosπlxcosπly2;m=211−a2;笔者试验,该公式是有效的,所以下文中我们将使用该公式.同时我们还可以考虑进行扫描方式的优化，上面方法的思路，是尽量使用新算出来的点的电势值.以下我们将使扫描尽量从边界出发，充分利用边界点电势初值的真实性。以下是就这个例子所做的尝试：当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值较为合适时，并不具有优势。由上文我们可以知道,这样的算法完全可以在计算机上有效实现,只要给定边界条件,就可以解出拉普拉斯方程或者泊松方程.我们现在还可以再使用这种算法解决和讨论一个实际的模型：相关问题：2.1模型引入:在《用计算机模拟闪电形成的尝试》--金秀儒，（以下称作文[1]）一文中，作者试图利用类似的方法解决闪电形成