第一类边界问题的有限差分法探讨文档格式.docx

1、三、 有限差分法公式的推导：把求解的区域划分成网格，把求解区域内连续的场分布用网络节点上的离散的数值解来代替。网格划分的充分细，才能够达到足够的精度。应用有限差分法计算静态场边值问题时，需要把微分方程用差分方程替代。用图形法解释如下:y卒有限差分网格划分如图将区域D划分为正方形网格，网格线的交点称为节点，两相邻平行网格线间的距离称为步距 h。然后，拉普拉斯方程离散化，对于任一点 0,有一阶偏导数：X=Xo（X。h,y。）- （X。- h,y。）2h而后，对于二阶偏导数:x（x。h/2,y。）- x（x。- h/2,y。x=x。h（p - 1 0 3h2对于丫轴同理:因此拉普拉斯方程的差分格式为

2、:+ P紧邻边界节点的拉普拉斯方程的差分格式为:亠亠p（1 p） q（1 q） 1 p 1 q其中p、q为小于1的正数；1、2为边界上的节点，其值为对应边界点处的值，是已知的。具体如图:4应用数值计算解释（泰勒公式展开法）：1点电位的泰勒公式展开为3点电位的泰勒公式展开为, 1 a跨平（訥方亦乔冰,当h很小时，忽略4阶以上的高次项，得同样可得将上面两式相加得口畀切丄护趴i pF = + I =- 在上式中代入 一，得肌=（F+約+血+的+如/4对于 1，即F=0的区域，得到二维拉普拉斯方程的有限差分形式旳=（幽+陷+伽十旳通过以上两种方法的推导，可得任意点的电位等于围绕它的四个点的 电位的平均

3、值。当用网格将区域划分后，对每一个网络点写出类似的 式子，就得到方程数与未知电位的网络点数相等的线性方程组。 已知的边界条件在离散化后成为边界上节点的已知电位值。四、差分方程组的解法 方法一：高斯一一赛德尔迭代法（简单迭代法） 其步骤是先对每一网格点设一初值。 然后按一个固定顺序（一般点的 顺序按“自然顺序”，即：从左到右，从下到上）如图:Ay56 7 89 9I_A之后，利用二维拉普拉斯方程的有限差分形式用围绕它的四个点的电 位的平均值作为它的新值，当所有的点计算完后，用它们的新值代替旧值，即完成了一次迭代。然后再进行下一次迭代，如此循环。如下式：i,r 1,2,:（k 1） _ 1 ： （

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

5、，把刚才计算得到的 临近点的新值代入，即在计算（i,j）点的电位时，把它左边的点（i-1,j）和 下面的点（i,j-1）的电位用刚才算过的新值代入，即：（迭代公式2） 第二，是引入“加速收敛因子”。上式中的a即为“加速收敛因子”， 且1a 2。特别关注的是逐次超松弛法收敛的快慢与 a有明显关 系。并且最佳a的取值随着条件的不同而不同，如何选择最佳 a，是 个复杂问题。在计算时可以尝试求取最佳 a值，以使计算快速。五、应用计算机仿真有限差分法解决具体问题 本次讨论我选择第四题为具体实例进行研究。题目如图所示，有一长方形的导体槽，a = 20，b = 15,设槽的长度为 无限长，槽上有一块与槽绝缘

6、的盖板，电位为 100V，其他板电 位为零，求槽内的电位分布。b=100 a x通过MATLAB进行仿真，运用有限差分法，源代码如下:u=zeros（15,20）;i=2:14;u（i,20）=100;for j=1:20for i=2:14u（i,j）=100/19*（j-1）;enda=in put（ please in put a（1a2）;a二 ）;m=i nput（ please in put m（1m）;m=for n=1:m;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

7、+f）;u（i,j）=b+a.*（g-b）;mesh（u）;通过对相关资料的查询及学习，我认为运用 MATLAB 编程就是 要将有限差分法的基本公式实现。 以此为出发点可以让解析和编写过 程相对简单化。首先，先要对电场的边界条件进行设定：使其满足题目的要求右侧边界电位为 100，其余三边为 0，由此 也可以判断出此题目为第一类边界条件的问题。同时，此处，也要特 别注意，因为在 MATLAB 程序中，编写一个矩阵的默认顺序是从左 到右，从上到下。另外，命名矩阵时是先编写行，后编写列，因此， 对于u（i,j）中，i代表的是第i行即对应纵坐标，j代表的是第j列即对 应横坐标，与坐标表示有一定的差别，

8、需要区分和注意。之后是相对关键的编写部分， 就是对电场内部的电位设定， 此处， 我通过分离变量法将u分解为u （x ）和u（y），分别求解。为让计算 简便，我设定其内电场为线性变化的，即 u（x） =kx+b， u（ y） =0 作为初值进行计算（应用分离变量法得到的结果是 u（x）二Bsinh mH/15 u（y）=Asin mH/15。可以通过此处在得出结果后进行对比和误差分析。） 又因为要在MATLAB进行运算，所以首先要转化为i和j的形式， 因此可得： u（j）=kj+b 和 u（i）=0 。从而得到方程组： u（1）=k+b=0;u（20）=20k+b=100。解得： u（j）=10

9、0/19*（j-1）, 即： u（i,j）=100/19*（j-1） 。由此得到了 内电场除去上下边界（衡为零）的电位方程，在 MATLAB中我运用了for 循环语句来实现其内电场的赋值，源代码如下： for j=1:将电场各点的电位赋初值完成后， 就需要带入公式进行计算了。 在之 确定。关于加速因子的最佳值确定问题，我通过学习和查阅相关资料得到了两种方法，如下：一是逐步搜索法；将a的取值区间（1,2 ）进行M等分，a分别取1+1/M,1+2/M,1+ （M-1） /M,通过上面提到的迭代公式2依次对同一个精度要求求出 迭代次数k的值，并从中选出“加速收敛因子” a的最佳值，具体步骤如下：步骤

10、1 给定等分数M和精度要求的值，令a的初始值为1；步骤2 令p=1,2,3, , M-1，重复步骤3-5 ;步骤 3 a p=1+p/M；步骤4 按照如下公式迭代ot.：（k 1） = .： （k） . _ : （k 1） . （k 1） . （k） . （k） _ 4 （k）i, j i,j i-1,j i,j -1 i 1, j i,j 1 4 i,j找出符合精度要求的迭代次数kp；步骤5 比较找出kp值最小的a p为最佳的加速收敛因子 a的值。二是黄金分割法：依据黄金分割法的思想，通过计算机自动选取最佳加速收敛因子a的近似值，具体步骤如下：步骤1 对（1,2 ）区间第一次0.618的分割

11、，区间分界a1=1,b1=2, 在（ a1, b）区间分割出黄金点point 1= a1+0.618（b 1- a 1）,进行迭代公 式2的迭代，求出迭代次数K值，如果迭代次数没有超出规定的发散常数，迭代结束，否则转为步骤 2步骤2 在（1,1.618 ）和（1.618,2 ）之间进行第二次的黄金 分割，找出分割点 point 2=a2+0.618（b2- a 2）,其中a?和b?是新分割区 间的左右边界。找出迭代次数最少的 a。以上两种方法是比较具有实践和研究价值进行的，在这里介绍和分享 给大家探讨。在我分析过后认为第二种较好，因为其可以通过计算机 自动选取到最佳的加速收敛因子，但是实施起来

12、比较复杂，要通过 MATLAB或其它语言进行编码实验得到，我初步试验过，没能成功， 就运用了第一种方法进行了计算，但是，我发现其运算过于复杂，笔 算正确的成功率比较低，而运用电脑编程要重复多次运算，虽然，计 算机运算快速且方便，但是，这样的源代码同样过于复杂，可以进一 步进行深入的研究。本次的题目解答中，由于是应用计算机计算，所以计算速度和计 算效率在不同的a取值中没有明显区别。因此，我就只是设定了一个 12之间的一个数字：1.8为“加速收敛因子” a的值了，进行了近 似计算，实验误差在可接受的范围。然后，是对于迭代次数的确定。当迭代满足停止条件：w （k*1） 申（k） wi,j i,j V

13、V时，对应的迭代次数即为最终的迭代次数， 这里的w为设定好的允许误差，为从简方式，我选择运用输入迭代次数，并通过多次输入对比 输出结果选取最佳迭代次数的方法来选取， 存在了一定的误差。但是, 相比于我运用上面的公式求解计算迭代次数， 进而得出的图形的误差 比要小许多， 我一直尝试进行改进， 但一直没能很好的解决误差过大 的问题：最主要的问题是边界值一直不满足我的设定值。 希望通过进 一步的研讨和大家的交流能过得到最优方式。通过以上的分析，我得出了进一步的编程代码：a=input（ please input a（1a=m=input（ please input m（1其中a和m分别为最佳“加速收

14、敛因子”和迭代次数，均为我们输入 数值计算的方式。先要预算出一些取值范围，然后通过分别输入，对 比输出结果进而来选取最优解。之后的源代码的编写思想就是通过迭代计算公式来实现对电场 内部各点电位值的计算，因此，在输入“加速收敛因子”和“迭代次 数”的数值后， 我通过 for 的循环语句实现了对电场内部的电位迭代 求解，并用图形形式表示出来。其源代码如下： for n=1: u（i,j）二b+a.*（g-b）;end end上述代码中，b,c,d,e,f 分别场内一点及其左右上下各点的电位值， 之后的g和u （i，j）则实现了逐次超松弛法的迭代，最后用 mesh 实现图形的输出。六、结果和误差的对

15、比分析由分离变量法得到的初步结果是u（x）二Bsinh mn /15和 u（y）=AsinmQ/15，从而应有 u=Csinh（mn/15）sinm（ n/15）的结果， 从其表达式的形式，我们可以简单推断出电场的电位在 x轴上的分布 是成虚指数函数的形式呈现，而在y轴上是成正弦函数呈现的，这和 通过有限差分法的计算机仿真的图形基本符合，如下图：因此，可以基本确定通过计算机仿真得到的结果的正确性。另外，可能存在较大误差的方面是跌代次数， 因为迭代次数是我 们人工输入的，所以判断其正确与否要通过输出的结果来判定、 来选 择，下面几幅图是当“加速收敛因子” a =1.8时，不同的迭代次数 m的对应

16、图形:m=1 时：Q if2Dm=10 时:m=20m=30m=32ifl.7m=40通过上面迭代次数 m取不同的值对应的不同图形分析可以看出，在 m=1,10,20时期计算结果都存在较大误差，图形与分离变量法得到的 结果相差很大，有明显的不规则部分，不成现光滑的曲线。而当m=30 后图形与分离变量法计算到的图形近似，而且图形规则，曲线光滑， 而且随着m值的增大图形的变化非常微小了，因此，我选择 m=32为迭代次数，所得到的结果误差比较小，结果相对准确。七、研究结论学习了运用有限差分法计算第一类边界条件的电位值， 通过计算机仿真进行了实际例题的求解，通过与分离变量法求解的对比，证明 可行性和有

17、效性，适用于解决第一类边界问题，而且更加方便和快捷。 同时，“加速收敛因子”最佳解问题比较复杂，不过通过计算机仿真 的过程中，不同数值“加速收敛因子”的加速效果不易体现出来。跌 代次数对计算结果又较大影响，迭代次数越多，结果越准确，误差越 小，但相应的计算也会更加复杂，因此只要迭代次数得到得结果满足 一定的允许误差范围即可。附录：1三类边界条件：第一类：已知电位在边界上的数值。第二类：已知电位的导数在边界上的数值。第三类：已知电位在部分边界上的数值和在其余边界上的导数的数值。2解析法包括：分离变量法，镜像法，复变函数法，保角变换法， 格林函数法。3数值法包括：数值积分法，有限差分法，有限元法，

18、矩量法。4解析法的优点：由有限个项构成闭合解或无穷级数，通常具有鲜 明的物理意义。缺点：解题范围窄小，对边界形状十分挑剔。数值法的优点：解题范围宽广，对边界的形状没有限制。 缺点：数据离散，物理意义深藏其中。心得体会：有限差分法的公式套用很简单， 记住结果很简单，但是要想深刻 理解其内容和灵活应用则比较困难， 需要一定时间的深刻学习和做题 应用，同时，这个方法确实是一个非常方便和有效地解决三类边界条 件问题的工具。MATLAB 在电磁场方面的应用相当广泛，其内拥有的很多函数 都可以模拟静电场中的实例， 就比如这次的研讨就应用了矩阵， 通过 有限差分法来模拟静电场的电位分布，相信将来还会有更多的有意 思，有研究价值计算和模拟电磁场的相关知识。 是学好电磁场必须要 牢固掌握的工具。