摘要本文针对解微分方程的有限元法,分别讨论了一维问.docVIP

摘要：本文针对解微分方程的有限元法，分别讨论了一维问题和二维问题，并给出了相应的程序实例，以便参考。由于面向对象程序设计的思想越发显得重要，所以本文的源程序也借用了这种思想，把有限元法看成是一个超类，而一维问题和二维问题都是在这个超类的基础上继承和发展起来的，所以程序结构清晰，便于阅读，而且代码冗余度低，扩充性好，例如扩展到三维问题原则上是可行的，只要理论支持。对于一维问题，本文分析了一个二阶椭圆方程，边值条件为第一和第二边值条件的混合问题，难度虽然不是很大，但对有限元解法的几个步骤给出了详尽的分析，例如变分问题，剖分问题，边值处理问题和误差分析问题，而且也给出了第三边值条件变分问题的具体表达式。对于二维问题，有限元法步骤同样，只是剖分问题和边值处理问题比较困难，所以本文以第一第二混合边值条件为例，对以上两个问题以专题的形式进行了讨论。由于二维问题选用了等边矩形作为剖分单元，而且原题又刚好是规则的正方形边界，所以剖分点编号相对容易，利用了一个二维数组记录剖分点，而在计算有限元方程组的时候又以行优先的原则把剖分点从左到右，从下到上进行了编号，至于边界处理问题，点的分类也相对容易，把第一边界上的点当作已知，减少了方程组未知数的个数，而第二边值条件也只是在某些方程组的右端多加了一个线积分。至于本实验跟别的组比较，优点就是采用的剖分元是二次矩形元而不是一次三角元，同等剖分点下，由于剖分单元本身固有的形状特性和剖分单元之间的对称性，计算积分次数也比三角元少，所以同等计算量下，精度要比其一次三角元高。总之，在故意提高难度和工作量的情况下，小组成员共同努力，相互鼓励支持使得本次实习完成得相当成功，既感谢老师给大家一个锻炼的机会，也感谢组内成员夜以继日的辛劳。 关键字：有限元法，面向对象，变分问题，剖分问题，边值处理，误差分析 实习题一（一维问题） 用线性元求下列边值问题的数值解： 为增加难度，考虑非齐次边值条件如下，只要令即得原问题的解。 ①问题分析： 所要求解的问题是两点边值问题的二阶常微分方程，而且是常系数非齐次方程，理论上来说总是可以通过先求对应的齐次方程的两个基解，再通过常数变易法得到非齐次方程的通解，最后代入边值条件便得到所要求的特解。而现在是通过变分法在有限维空间得到其数值解，即有限元分析法，从而可以避免像求理论解那样为求表达式而进行复杂的积分运算，而且有些积分是暂时无法用解析式来表达的，最终还是得通过数值方法来解决。 用解常系数二阶常微分方程方法得到上面方程的理论解的解析式为： ②假设： 这里总是假设计算机的计算是精确的，即无舍入误差，不需考虑绝对稳定和相对稳定性 ③建立模型： 对一般的两点边值问题： 利用虚功原理得到其变分形式为 (*) 其中， 即求满足边值条件的，对任给的满足边值条件的v使变分方程(*)恒成立。 ④求解： 用Glalerkin法，难度增加题中所给条件对应模型中 相应的变分形式为 2，对区间（0，1）进行等长剖分，假设其剖分点个数包括两端点为n+1个，即剖分点为，其步长h是随n而变的。 为每个剖分点选取线性插值函数作为基函数，即 4，由基函数计算刚度矩阵元素和非齐次项系数，并根据边界条件适当处理某些系数（具体实现算法在源程序代码中给出）； 5，解有限元方程组 ⑤上机实现：请对照源程序阅读 为便于讨论，以下分析如不特别说明，则取N＝10， 剖分 按照文件输入的N，剖分点个数为N＋1，且计算公式 ： height=high-low; h=steplen=height/N; t[i]=low+i*steplen;（其中low＝0,high=1） 计算有限元方程组系数 for(j=1;jN;j++){ stiffness[j-1][j]=GL(f1,0,1,4); stiffness[j][j]=GL(f2,0,1,4)+GL(f3,0,1,4); stiffness[j+1][j]=GL(f1,0,1,4); } stiffness[j][j]就是所要计算的刚度矩阵系数，形如GL(f2,0,1,4)的函数为高斯勒让德积分公式，第一个参数对应基函数，第二三个参数为积分上下限，第四个参数为（点数－1），记录积分公式所取点数，这里都是用4，即五点公式。请注意当j＝N的时候要单独考虑： stiffness[N-1][N]=GL(f1,0,1,4); stiffness[N][N]=GL(f2,0,1,4); 因为最后一个基函数只跨一个剖分区间，而前面的基函数除了外都占了两个剖分区间，这个在下面的三对角系数刚度矩阵中可以看到，表现为最后一个对角矩阵元素的值只有前面对角元素值的1/2