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

1、 第一类边界问题的有限差分法探讨 摘要：本次重点是对于第一类边界问题的两种不同方法的比照研讨， 通过计算机仿真有限差分法和计算别离变量法对同一问题的求解， 对 结果进行比照，能够发现有限差分法更加快捷简便， 只要迭代次数足 够多就能使误差趋于零。而别离变量法那么是准确的计算出结果， 只是 运算相对复杂。 关键字：有限差分法，别离变量法，加速收敛因子，迭代次数，边界 条件。 引言：在给定的三类边界条件下求解标量位或矢量位的泊松方程或 拉普拉斯方程的解一般的理论依据是唯一性定理和得加原理， 由此而 得出的解题方法有很多。主要分为两大类：一是解析法如别离变量 法，镜像法等，二是数值法如有限差分法，有

2、限元法 等。这两 种方法各有优点和缺乏，相比拟而言在许多实际问题中由于边界条 件过于复杂而无法求得解析解。这就需要借助于数值法来求电磁场的 数值解。有限差分法便是一种比拟容易的数值解法。 本次研讨就以第 一类边界问题进行为例来分析研究有限差分法。 一、有限差分法的定义： 微分方程和积分微分方程数值解的方法为有限差分法。 根本 思想是把连续的定解区域用有限个离散点构成的网格来代替， 这 些离散点称作网格的节点；把连续定解区域上的连续变量的函数 用在网格上定义的离散变量函数来近似；把原方程和定解条件中 的微商用差商来近似， 积分用积分和来近似， 于是原微分方程和 定解条件就近似地代之以代数方程组，

3、即有限差分方程组 ，解 此方程组就可以得到原问题在离散点上的近似解。然后再利用插 值方法便可以从离散解得到定解问题在整个区域上的近似解。 二、 有限差分法解题的根本步骤： (1) 、区域离散化，即把所给偏微分方程的求解区域细分成由有 限个格点组成的网格； (2) 、近似替代，即采用有限差分公式替代每一个格点的导数； (3) 、逼近求解。换而言之，这一过程可以看作是用一个插值多 项式及其微分来代替偏微分方程的解的过程。 三、 有限差分法公式的推导： 把求解的区域划分成网格，把求解区域内连续的场分布用 网络节点上的离散的数值解来代替。网格划分的充分细，才能 够到达足够的精度。应用有限差分法计算静态

4、场边值问题时， 需要把微分方程用差分方程替代。 用图形法解释如下： 有限差分网格划分 在由边界L界定的二维区域 D内，电位函数满足拉普拉斯 方程且给定第一边界条件，贝 2 ：2 - - =0 2 2 x y * IL= f(s) 如图将区域D划分为正方形网格，网格线的交点称为节点， 两相邻平行网格线间的距离称为步距 h。然后，拉普拉斯方程 离散化，对于任一点 0,有一阶偏导数： (xo h,y。)- (x- h,yo)=, ：Xx=xo 2h x 而后，对于二阶偏导数: x(xo h/2,y。)- x(x。- h/2,yo) h 中-+中 1 0 3(在区域D内) x=xo p 2 因此拉普拉

5、斯方程的差分格式为: +平 +中 +中 -A? 1 2 3 4 4 0 紧邻边界节点的拉普拉斯方程的差分格式为: 乳 八 八 中 (i i _ 1 _ + _2 _ + + _ 4 _ + _ =0 p(1 + p) q(1 + q) 1 + p 1 + q p qj 0 其中p、q为小于1的正数；1、2为边界上的节点，其值为对应边界 点处的值，是的。具体如图: 2 qh 3 1 h ph h 4 应用数值计算解释(泰勒公式展开法)： 1点电位的泰勒公式展开为 对于Y轴同理: P 4 紧邻边界节点的网格划分 ,耶 1 ,3 0、. 2 1诲伊3 + (+Tr)洲 * 海 2! ax2 3!=

6、3点电位的泰勒公式展开为 海 1 8 2 1 .8申皿 /。一京小融蒙。与京。s + =? + f* L 灼期饥a? ，当h很小时，忽略4阶以上的高次 项，得 并0炉注孔-2饥 同样可得 芋冲=色+知一 2险 将上面两式相加得 a% 7 + 570 =代+代+灼+知一机0 4W上 在上式中代入 脆朗 % ,得 所=用护+啊+ % +的+R/4 对于。=,即F=0的区域，得到二维拉普拉斯方程的有限差分形式f 寸2& +， 仞广的+啊+的+也苗 通过以上两种方法的推导，可得任意点的电位等于围绕它的四个点的 电位的平均值。当用网格将区域划分后，对每一个网络点写出类似的 式子，就得到方程数与未

7、知电位的网络点数相等的线性方程组。 的边界条件在离散化后成为边界上节点的电位值。 四、差分方程组的解法 方法一：高斯一一赛德尔迭代法简单迭代法 其步骤是先对每一网格点设一初值。 然后按一个固定顺序一般点的 顺序按“白然顺序，即：从左到右，从下到上如图: y 5 _ 4 6 _ J 7 _ j 8 _ 9 |_24 9 V |44 1_ . o 1 X 之后，利用二维拉普拉斯方程的有限差分形式用围绕它的四个点的电 式: (k 1) . (k) . (k) i,j -1 i 1,j i,j 1 ，j 位的平均值作为它的新值, 当所有的点计算完后，用它们的新值代替 旧值，即完成了一次迭代。 然后再进

8、行下一次迭代，如此循环。如下 ，(k 1) _ 1 .= (k 1). i,j 4 i-1,j 迭代公式1 其式中的上角标k表示k次近似值，下脚标i, j表示节点所在位 置，即第i行第j列的交点。其中要特别注意：在迭代过程中遇到边 界点式，需将边界条件 中.=f.带入。 循环迭代时当所有内节点满足以下条件时停止迭代： (k 1) i,j 其中，W是预定的最大允许误差。 方法二：逐次超松弛法 简单迭代法在解决问题时收敛速度比拟慢， 实用价值不大。实际中常 采用逐次超松弛法又称高斯赛德尔迭代法变形，相比之下它 有两点重大的改良，第一是计算每一网格点时，把刚刚计算得到的 临近点的新值代入，即在计算i

9、,j点的电位时，把它左边的点i-1,j和 下面的点i,j-1的电位用刚刚算过的新值代入，即： (k 1) _ (k).二：(k 1) . ： (k 1) . ： (k) . (k) (k) i,j i,j 4 i-1,j i,j-1 i 1,j i,j 1 4 i,j 迭代公式2 第二，是引入“加速收敛因子。上式中的a即为“加速收敛因子 且1 a 2。特别关注的是逐次超松弛法收敛的快慢与 a有明显关 系。并且最正确a的取值随着条件的不同而不同，如何选择最正确 a ,是 个复杂问题。在计算时可以尝试求取最正确 a值，以使计算快速。 五、应用计算机仿真有限差分法解决具体问题 本次讨论我选择第四题为

10、具体实例进行研究。 题目： (k) i,j 如下图，有一长方形的导体槽，a = 20, b = 15,设槽的长度为 无限长，槽上有一块与槽绝缘的盖板，电位为 100V,其他板电 位为零，求槽内的电位分布。 b =100 - - a x 通过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); end end a=input( please input a(1a2);a= ); m=input( please input m(1m);m=

11、); 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 end end mesh(u); 通过对相关资料的查询及学习，我认为运用 MATLAB编程就是 要将有限差分法的根本公式实现。以此为出发点可以让解析和编写过 程相对简单化。 首先，先要对电场的边界条件进行设定： u=zeros(15,20); i=2:14; u(i,20)=100; 使其满足题目的要求右侧边界电位为 10

12、0,其余三边为0,由此 也可以判断出此题目为第一类边界条件的问题。同时，此处，也要特 别注意，因为在MATLAB程序中，编写一个矩阵的默认顺序是从左 到右，从上到下。另外，命名矩阵时是先编写行，后编写列，因此， 对于u(i,j)中，i代表的是第i行即对应纵坐标，j代表的是第j列即对 应横坐标，与坐标表示有一定的差异，需要区分和注意。 之后是相对关键的编写局部，就是对电场内部的电位设定，此处， 我通过别离变量法将u分解为u (x)和u(y),分别求解。为让计算 简便，我设定其内电场为线性变化的，即 u (x) =kx+b, u (y) =0 作为初值进行计算(应用别离变量法得到的结果是 u(x)

13、=Bsinh nf!/15 u(y)=Asin nH/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); e

14、nd end 将电场各点的电位赋初值完成后，就需要带入公式进行计算了在之 确定。关于加速因子的最正确值确定问题， 我通过学习和查阅相关资料 得到了两种方法，如下： 一是逐步搜索法； 将a的取值区间(1,2 )进行M等分，a分别取1+1/M,1+2/M, 1+ (M-1) /M,通过上面提到的迭代公式2依次对同一个精度要求求出 迭代次数k的值，并从中选出“加速收敛因子 a的最正确值，具体 步骤如下： 步骤1 给定等分数M和精度要求的值，令a的初始值为1； 步骤2 令p=1,2,3, , M-1，重复步骤3-5; 步骤3 以p=1+p/M; 步骤4 按照如下公式迭代 . 、 a . . 八、 八、

15、 、 :(k 1).： (k) . _ .： (k 1) . (k 1) . (k) . . (k) . (k) i,j i,j 4 i-1,j i,j-1 i 1,j i,j 1 4 i,j 找出符合精度要求&的迭代次数kp； 步骤5 比拟找出kp值最小的以p为最正确的加速收敛因子 a的值。 p 二是黄金分割法： 依据黄金分割法的思想，通过计算机白动选取最正确加速收敛因子 a的近似值，具体步骤如下： 步骤1 对(1,2 )区间第一次0.618的分割，区间分界a1=1,b1=2, 在(a, b)区间分割出黄金点point 1= a+0.618(b 1- a 1),进行迭代公 前先要对“

16、加速收敛因子 a的最正确值和预定的最大允许误差进行 式2的迭代，求出迭代次数K值，如果迭代次数没有超出规定的发散 常数，迭代结束，否那么转为步骤 2。 步骤2 在1,1.618 和1.618,2 之间进行第二次的黄金 分割，找出分割点 point 2=a2+0.618b2- a 2,其中a2和b2是新分割区 间的左右边界。找出迭代次数最少的 a。 以上两种方法是比拟具有实践和研究价值进行的，在这里介绍和分享 给大家探讨。在我分析过后认为第二种较好，因为其可以通过计算机 白动选取到最正确的加速收敛因子，但是实施起来比拟复杂，要通过 MATLAB其它语言进行编码实验得到，我初步试验过，没能成功，

17、就运用了第一种方法进行了计算，但是，我发现其运算过于复杂，笔 算正确的成功率比拟低，而运用电脑编程要重复屡次运算，虽然，计 算机运算快速且方便，但是，这样的源代码同样过于复杂，可以进一 步进行深入的研究。 本次的题目解答中，由于是应用计算机计算，所以计算速度和计 算效率在不同的a取值中没有明显区别。因此，我就只是设定了一个 12之间的一个数字：1.8为“加速收敛因子 a的值了，进行了近 似计算，实验误差在可接受的范围。 然后，是对于迭代次数确实定。当迭代满足停止条件： (k 1) _ (k) . . . i,j i,j 时，对应的迭代次数即为最终的迭代次数， 这里的W为设定好的允许 误差，为从

18、简方式，我选择运用输入迭代次数，并通过屡次输入比照 输出结果选取最正确迭代次数的方法来选取， 存在了一定的误差。但是, 相比于我运用上面的公式求解计算迭代次数， 进而得出的图形的误差 比要小许多，我一直尝试进行改良，但一直没能很好的解决误差过大 的问题：最主要的问题是边界值一直不满足我的设定值。 希望通过进 一步的研讨和大家的交流能过得到最优方式。 通过以上的分析，我得出了进一步的编程代码： a=input( please input a(1a2);a= ); m=input( please input m(1m);m= ); 其中a和m分别为最正确“加速收敛因子和迭代次数，均为我们输入 数值

19、计算的方式。先要预算出一些取值范围，然后通过分别输入，对 比输出结果进而来选取最优解。 之后的源代码的编写思想就是通过迭代计算公式来实现对电场 内部各点电位值的计算，因此，在输入“加速收敛因子和“迭代次 数的数值后，我通过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 end end me

20、sh(u); 上述代码中，b,c,d,e,f 分别场内一点及其左右上下各点的电位值， 之后的g和u (i , j)那么实现了逐次超松弛法的迭代，最后用 mesh 实现图形的输出。 六、结果和误差的比照分析 由别离变量法得到的初步结果是u(x)=Bsinh mEt /15和 u(y)=Asinnf!/15 ,从而应有 u=Csinh(mH/15)sinm( n/15)的结果， 从其表达式的形式，我们可以简单推断出电场的电位在 x轴上的分布 是成虚指数函数的形式呈现，而在y轴上是成正弦函数呈现的，这和 通过有限差分法的计算机仿真的图形根本符合，如下列图： 因此，可以根本确定通过计算机仿真得到的结果

21、的正确性。 另外，可能存在较大误差的方面是跌代次数， 因为迭代次数是我 们人工输入的，所以判断其正确与否要通过输出的结果来判定、 来选 择，下面几幅图是当“加速收敛因子 a =1.8时，不同的迭代次数 m的对应图形: m=1 时: m=10 时: m=20 皿 m=30 m=32 却417 m=35 m=40 通过上面迭代次数 m取不同的值对应的不同图形分析可以看出，在 m=1,10,20时期计算结果都存在较大误差，图形与别离变量法得到的 结果相差很大，有明显的不规那么局部，不成现光滑的曲线。而当m=30 后图形与别离变量法计算到的图形近似，而且图形规那么，曲线光滑， 而且随着m值的增大图形的变化非常微小了，因此，我选择 m=32 为迭代次数，所得到的结果误差比拟小，结果相对准确。 七、研究结论 学习了运