第一类边界问题的有限差分法探讨.doc

第一类边界问题的有限差分法探讨摘要：本次重点是对于第一类边界问题的两种不同方法的对比研讨，通过计算机仿真有限差分法和计算分离变量法对同一问题的求解，对结果进行对比，能够发现有限差分法更加快捷简便，只要迭代次数足够多就能使误差趋于零。而分离变量法则是准确的计算出结果，只是运算相对复杂。关键字:有限差分法，分离变量法，加速收敛因子，迭代次数，边界条件。引言：在给定的三类边界条件下求解标量位或矢量位的泊松方程或拉普拉斯方程的解一般的理论依据是唯一性定理和得加原理，由此而得出的解题方法有很多。主要分为两大类：一是解析法（如分离变量法，镜像法等），二是数值法（如有限差分法，有限元法等）。这两种方法各有优点和不足，相比较而言在许多实际问题中由于边界条件过于复杂而无法求得解析解。这就需要借助于数值法来求电磁场的数值解。有限差分法便是一种比较容易的数值解法。本次研讨就以第一类边界问题进行为例来分析研究有限差分法。一、 有限差分法的定义：微分方程和积分微分方程数值解的方法为有限差分法。基本思想是把连续的定解区域用有限个离散点构成的网格来代替， 这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似， 积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组 ， 解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。二、 有限差分法解题的基本步骤：（1）、区域离散化，即把所给偏微分方程的求解区域细分成由有限个格点组成的网格； （2）、近似替代，即采用有限差分公式替代每一个格点的导数； （3）、逼近求解。换而言之，这一过程可以看作是用一个插值多项式及其微分来代替偏微分方程的解的过程。三、 有限差分法公式的推导：把求解的区域划分成网格，把求解区域内连续的场分布用网络节点上的离散的数值解来代替。网格划分的充分细，才能够达到足够的精度。应用有限差分法计算静态场边值问题时，需要把微分方程用差分方程替代。用图形法解释如下：有限差分网格划分在由边界L界定的二维区域D内，电位函数满足拉普拉斯方程且给定第一边界条件，则：如图将区域D划分为正方形网格，网格线的交点称为节点，两相邻平行网格线间的距离称为步距 h。然后，拉普拉斯方程离散化，对于任一点0，有一阶偏导数：而后，对于二阶偏导数：对于Y轴同理：因此拉普拉斯方程的差分格式为：紧邻边界节点的拉普拉斯方程的差分格式为：其中p、q为小于1的正数；1、2为边界上的节点，其值为对应边界点处的值，是已知的。具体如图：紧邻边界节点的网格划分应用数值计算解释（泰勒公式展开法）：1点电位的泰勒公式展开为3点电位的泰勒公式展开为，当h很小时，忽略4阶以上的高次项，得 同样可得将上面两式相加得在上式中代入，得对于，即F=0的区域，得到二维拉普拉斯方程的有限差分形式通过以上两种方法的推导，可得任意点的电位等于围绕它的四个点的电位的平均值。当用网格将区域划分后，对每一个网络点写出类似的式子，就得到方程数与未知电位的网络点数相等的线性方程组。已知的边界条件在离散化后成为边界上节点的已知电位值。四、 差分方程组的解法方法一：高斯赛德尔迭代法（简单迭代法）其步骤是先对每一网格点设一初值。然后按一个固定顺序（一般点的顺序按“自然顺序”，即：从左到右，从下到上）如图：之后，利用二维拉普拉斯方程的有限差分形式用围绕它的四个点的电位的平均值作为它的新值，当所有的点计算完后，用它们的新值代替旧值，即完成了一次迭代。然后再进行下一次迭代，如此循环。如下式：（迭代公式1）其式中的上角标（k）表示k次近似值，下脚标i，j表示节点所在位置，即第i行第j列的交点。其中要特别注意：在迭代过程中遇到边界点式，需将边界条件 带入。循环迭代时当所有内节点满足以下条件时停止迭代：其中，W是预定的最大允许误差。方法二：逐次超松弛法 简单迭代法在解决问题时收敛速度比较慢，实用价值不大。实际中常采用逐次超松弛法（又称高斯赛德尔迭代法变形），相比之下它有两点重大的改进 ，第一是计算每一网格点时，把刚才计算得到的临近点的新值代入，即在计算(i,j)点的电位时，把它左边的点(i-1,j)和下面的点(i,j-1)的电位用刚才算过的新值代入，即:（迭代公式2）第二，是引入“加速收敛因子”。上式中的即为“加速收敛因子”，且1 2。特别关注的是逐次超松弛法收敛的快慢与有明显关系。并且最佳的取值随着条件的不同而不同，如何选择最佳，是个复杂问题。在计算时可以尝试求取最佳值，以使计算快速。五、 应用计算机仿真有限差分法解决具体问题本次讨论我选择第四题为具体实例进行研究。题目：如图所示，有一长方形的导体槽，a = 20，b = 15，设槽的长度为无限长，槽上有一块与槽绝缘的盖板，电位为100V，其他板电位为零，求槽内的电位分布。 通过MATLAB进行仿真，运用有限差分法，源代码如下：u=zeros(15,20);i=2:14;u(i,20)=100;for j=1:20 for i=2:14 u(i,j)=100/19*(j-1); endenda=input(please input a(1a2);a=);m=input(please input m(1m);m=);for n=1:m; for i=2:14; for j=2:19; b=u(i,j); c=u(i,j+1); d=u(i+1,j); e=u(i-1,j); f=u(i,j-1); g=0.25.*(c+d+e+f); u(i,j)=b+a.*(g-b); end endendmesh(u);通过对相关资料的查询及学习，我认为运用MATLAB编程就是要将有限差分法的基本公式实现。以此为出发点可以让解析和编写过程相对简单化。首先，先要对电场的边界条件进行设定：u=zeros(15,20);i=2:14;u(i,20)=100;使其满足题目的要求右侧边界电位为100，其余三边为0，由此也可以判断出此题目为第一类边界条件的问题。同时，此处，也要特别注意，因为在MATLAB程序中，编写一个矩阵的默认顺序是从左到右，从上到下。另外，命名矩阵时是先编写行，后编写列，因此，对于u(i,j)中，i代表的是第i行即对应纵坐标，j代表的是第j列即对应横坐标，与坐标表示有一定的差别，需要区分和注意。之后是相对关键的编写部分，就是对电场内部的电位设定，此处，我通过分离变量法将u分解为u（x）和u(y)，分别求解。为让计算简便，我设定其内电场为线性变化的，即u（x）=kx+b，u（y）=0作为初值进行计算（应用分离变量法得到的结果是u(x)=Bsinh m/15u(y)=Asinm/15。可以通过此处在得出结果后进行对比和误差分析。)又因为要在MATLAB中进行运算，所以首先要转化为i和j的形式，因此可得：u(j)=kj+b和u(i)=0。从而得到方程组：u(1)=k+b=0;u(20)=20k+b=100。解得：u(j)=100/19*(j-1),即：u(i,j)=100/19*(j-1)。由此得到了内电场除去上下边界（衡为零）的电位方程，在MATLAB中我运用了for循环语句来实现其内电场的赋值，源代码如下：for j=1:20 for i=2:14 u(i,j)=100/19*(j-1); endend将电场各点的电位赋初值完成后，就需要带入公式进行计算了。在之前先要对“加速收敛因子” 的最佳值和预定的最大允许误差进行确定。关于加速因子的最佳值确定问题，我通过学习和查阅相关资料得到了两种方法，如下：一是逐步搜索法；将的取值区间（1,2）进行M等分，分别取1+1/M,1+2/M1+（M-1）/M,通过上面提到的迭代公式2依次对同一个精度要求求出迭代次数k的值，并从中选出“加速收敛因子” 的最佳值，具体步骤如下：步骤1 给定等分数M和精度要求的值，令的初始值为1；步骤2 令p=1,2,3,M-1，重复步骤3-5；步骤3 p=1+p/M;步骤4 按照如下公式迭代找出符合精度要求的迭代次数kp;步骤5 比较找出kp值最小的p为最佳的加速收敛因子的值。二是黄金分割法：依据黄金分割法的思想，通过计算机自动选取最佳加速收敛因子的近似值，具体步骤如下：步骤1 对（1,2）区间第一次0.618的分割，区间分界a1=1,b1=2,在（a1，b1）区间分割出黄金点point1= a1+0.618(b1- a1),进行迭代公式2的迭代，求出迭代次数K值，如果迭代次数没有超出规定的发散常数，迭代结束，否则转为步骤2。步骤2 在（1,1.618）和（1.618,2）之间进行第二次的黄金分割，找出分割点point2=a2+0.618(b2- a2),其中a2和b2是新分割区间的左右边界。找出迭代次数最少的。以上两种方法是比较具有实践和研究价值进行的，在这里介绍和分享给大家探讨。在我分析过后认为第二种较好，因为其可以通过计算机自动选取到最佳的加速收敛因子，但是实施起来比较复杂，要通过MATLAB或其它语言进行编码实验得到，我初步试验过，没能成功，就运用了第一种方法进行了计算，但是，我发现其运算过于复杂，笔算正确的成功率比较低，而运用电脑编程要重复多次运算，虽然，计算机运算快速且方便，但是，这样的源代码同样过于复杂，可以进一步进行深入的研究。本次的题目解答中，由于是应用计算机计算，所以计算速度和计算效率在不同的取值中没有明显区别。因此，我就只是设定了一个12之间的一个数字：1.8为“加速收敛因子” 的值了，进行了近似计算，实验误差在可接受的范围。然后，是对于迭代次数的确定。当迭代满足停止条件：时，对应的迭代次数即为最终的迭代次数，这里的W为设定好的允许误差，为从简方式，我选择运用输入迭代次数，并通过多次输入对比输出结果选取最佳迭代次数的方法来选取，存在了一定的误差。但是，相比于我运用上面的公式求解计算迭代次数，进而得出的图形的误差比要小许多，我一直尝试进行改进，但一直没能很好的解决误差过大的问题：最主要的问题是边界值一直不满足我的设定值。希望通过进一步的研讨和大家的交流能过得到最优方式。通过以上的分析，我得出了进一步的编程代码：a=input(please input a(1a2);a=);m=input(please input m(1m);m=);其中a和m分别为最佳“加速收敛因子”和迭代次数，均为我们输入数值计算的方式。先要预算出一些取值范围，然后通过分别输入，对比输出结果进而来选取最优解。之后的源代码的编写思想就是通过迭代计算公式来实现对电场内部各点电位值的计算，因此，在输入“加速收敛因子”和“迭代次数”的数值后，我通过for的循环语句实现了对电场内部的电位迭代求解，并用图形形式表示出来。其源代码如下：for n=1:m; for i=2:14; for j=2:19; b=u(i,j); c=u(i,j+1); d=u(i+1,j); e=u(i-1,j); f=u(i,j-1); g=0.25.*(c+d+e+f); u(i,j)=b+a.*(g-b); end endendmesh(u);上述代码中，b,c,d,e,f分别场内一点及其左右上下各点的电位值，之后的g和u（i，j）则实现了逐次超松弛法的迭代，最后用mesh实现图形的输出。六、 结果和误差的对比分析由分离变量法得到的初步结果是u(x)=Bsinh m/15和u(y)=Asinm/15，从而应有u=Csinh(m/15)sinm(/15)的结果，从其表达式的形式，我们可以简单推断出电场的电位在x轴上的分布是成虚指数函数的形式呈现，而在y轴上是成正弦函数呈现的，这和通过有限差分法的计算机仿真的图形基本符合，如下图：因此，可以基本确定通过计算机仿真得到的结果的正确性。另外，可能存在较大误差的方面是跌代次数，因为迭代次数是我们人工输入的，所以判断其正确与否要通过输出的结果来判定、来选择，下面几幅图是当“加速收敛因子” =1.8时，不同的迭代次数m的对应图形：m=1时：m=10时：m=20m=30m=32m=35m=40通过上面迭代次数m取不同的值对应的不同图形分析可以看出，在m=1,10,20时期计算结果都存在较大误差，图形与分离变量法得到的结果相差很大，有明显的不规则部分，不成现光滑的曲线。而当m=30后图形与分离变量法计算到的图形近似，而且图形规则，曲线光滑，而且随着m值的增大图形的变化非常微小了，因此，我选择m=32为迭代次数，所得到的结果误差比较小，结果相对准确。七、 研究结论学习了运用有限差分法计算第一类边界条件的电位值，通过计算机仿真进行了实际例题的求解，通过与分离变量法求解的对比，证明可行性和有效性，适用于解决第一类边界问题，而且更加方便和快捷。同时，“加速收敛因子”最佳解问题比较复杂，不过通过计算机仿真的过程中，不同数值“加速收敛因子”的加速效果不易体现出来。跌代次数对计算结果又较大影响，迭代次数越多，结果越准确，误差越小，但相应的计算也会更加复杂，因此只要迭代次数得到得结果满足一定的允许误差范围即可。附录： 三类边界条件：第一类：已知电位在边界上的数值。第二类：已知电位的导数在边界上的数值。第三类：已知电位在一部分边界上的数值和在其余边界上的导数的数值。 解析法包括：分离变量法，镜像法，复变函数法，保角变换法，格林函数法。 数值法包括：数值积分法，有限差分法，有限元法，矩量法。 解析法的优点：由有限个项构成闭合解或无穷级数，通常具有鲜