第25卷第7期 重庆理工大学学报（自然科学） 2011年7月 Vol.25 No.7 Journal of Chongqing University of Technology(Natural Science) Jul.2011 复连通二维Stokes问题的Galerkin边界元解法 王小军，刘晓宇，杜亚楠 (重庆大学数学与统计学院，重庆400030) 摘 要：介绍了复连通区域二维Stokes问题的Galerkin边界元解法的思路与程序设计方 法。利用Fortran计算所求，点的流速，用Matlab画出相应的流线图。在计算系数时，对第1重积 分中存在奇异性的单元采用解析积分，不存在奇异性的单元采用数值积分，第2重不存在奇异 性，故采用数值积分。用若干数值算例模拟了复连通区域上不可压缩黏性流体的绕流。 关键词：Stokes问题；边界元；复连通；程序设计 中图分类号：0241.82 文献标识码：A 文章编号：1674-8425(2011)07-0117-04 BEM for Solving Stokes Problem in Complex Connected Domain WANG Xiao-jun,LIU Xiao-yu,DU Ya-nan (College of Mathematics and Statistics,Chongqing University,Chongqing 400030,China) Abstract:In this paper,we deduce the boundary integral equation of 2-Stokes problem in a com- plex connected region.Then we solve the boundary variational equation by Galerkin boundary element method.Velocity of flow is calculated by Fortran programme while streamlined diagram is drawn by Matlab programme.During the course of computing coefficient,numerical integral is used when there is no singularity while analytic integral is used when there is singularity.In this paper,several numer- ical tests simulate the viscous incompressible flow in a complex connected region. Key words:stokes problem;boundary element method;complex connected domain;programming 边界元方法是求解不可压缩粘性流体Stokes 计等技术性问题进行探讨。 流动问题数值解的理想方法。对于二维问题，在 根据Green公式和Stokes算子的基本解可以 区域的边界是封闭曲线的情况下，如无界流体的 推导出对应于复连通闭曲线Stokes问题的边界积 绕流问题或有界固壁容器中流体的流动问题等 分方程，得出基于单层位势的间接边界积分方程， Stokes问题，用边界元方法求解已有丰富的研究成 这是一个第1类的Fredholm向量积分方程。对与 果，然而多侧重于理论分析。本文拟对复连通二 之等价的边界变分方程采用Galerkin边界元求 维Stokes问题的Galerkin边界元解法中的程序设 解，得出单层位势的向量密度，进而得出流场中任 收稿日期：2010-11-23 基金项目：国家十一五科技支撑计划(2006BAU02A09) 作者简介：王小军(1984一），男，山东青州人，硕士研究生，主要从事边界元及其在工程计算中的应用研究。 ?1994-2017 China Academic Journal Electronic Publishing House.All rights reserved.http://www.cnki.net
第 25 卷 第 7 期 Vol． 25 No． 7 重庆理工大学学报( 自然科学) Journal of Chongqing University of Technology( Natural Science) 2011 年 7 月 Jul． 2011 收稿日期:2010 － 11 － 23 基金项目:国家十一五科技支撑计划( 2006BAJ02A09) 作者简介:王小军( 1984—) ，男，山东青州人，硕士研究生，主要从事边界元及其在工程计算中的应用研究。 复连通二维 Stokes 问题的 Galerkin 边界元解法 王小军，刘晓宇，杜亚楠 ( 重庆大学 数学与统计学院，重庆 400030) 摘 要:介绍了复连通区域二维 Stokes 问题的 Galerkin 边界元解法的思路与程序设计方 法。利用 Fortran 计算所求点的流速，用 Matlab 画出相应的流线图。在计算系数时，对第 1 重积 分中存在奇异性的单元采用解析积分，不存在奇异性的单元采用数值积分，第 2 重不存在奇异 性，故采用数值积分。用若干数值算例模拟了复连通区域上不可压缩黏性流体的绕流。 关 键 词:Stokes 问题; 边界元; 复连通; 程序设计 中图分类号:O241． 82 文献标识码:A 文章编号:1674 － 8425( 2011) 07 － 0117 － 04 BEM for Solving Stokes Problem in Complex Connected Domain WANG Xiao-jun，LIU Xiao-yu，DU Ya-nan ( College of Mathematics and Statistics，Chongqing University，Chongqing 400030，China) Abstract: In this paper，we deduce the boundary integral equation of 2-D Stokes problem in a complex connected region． Then we solve the boundary variational equation by Galerkin boundary element method． Velocity of flow is calculated by Fortran programme while streamlined diagram is drawn by Matlab programme． During the course of computing coefficient，numerical integral is used when there is no singularity while analytic integral is used when there is singularity． In this paper，several numerical tests simulate the viscous incompressible flow in a complex connected region． Key words: stokes problem; boundary element method; complex connected domain; programming 边界元方法是求解不可压缩粘性流体 Stokes 流动问题数值解的理想方法。对于二维问题，在 区域的边界是封闭曲线的情况下，如无界流体的 绕流问题或有界固壁容器中流体的流动问题等 Stokes 问题，用边界元方法求解已有丰富的研究成 果，然而多侧重于理论分析。本文拟对复连通二 维 Stokes 问题的 Galerkin 边界元解法中的程序设 计等技术性问题进行探讨。 根据 Green 公式和 Stokes 算子的基本解可以 推导出对应于复连通闭曲线 Stokes 问题的边界积 分方程，得出基于单层位势的间接边界积分方程， 这是一个第 1 类的 Fredholm 向量积分方程。对与 之等价的边界变分方程采用 Galerkin 边界元求 解，得出单层位势的向量密度，进而得出流场中任
118 重庆理工大学学报 意点的流速值。本文主要讨论其中用利用Fortran 的关系，由于“在穿越边界时的连续性，从式(2) 计算所求点的流速与用Matlab画出相应的流线图 得到联系己知函数。和中间变量t的积分方 的主要过程。 程组： 1复连通区域Stokes问题的边界积分方程 uot(y）= k=1,2,y∈T (3) 设Γ是R中有限区域2的封闭曲线，2是 对于方程(3)，方程两边同乘以1(y),然后在T上 T形成的封闭区域中的一条封闭曲线，形成的 积分，便得到变分方程组： 区域为2,2=2-2。设T=+2,考虑如 下的问题： 2iypu.y,ed山，= r-v△u+gradp=0,在2中 (4) div u =0, 在2中 (1) 含s (urr=uo, 在T+2上 2.2离散变分方程 为了数值求解变分方程(4)，边界T必须离 其中：u=(山1，山2)是流速；p是压强；u。是闭曲线边 散为一系列单元，从而把边界积分方程转化为线 界T以及内部闭曲线边界2上的己已知函数。 性代数方程组进行求解。设把边界T离散为 个单元，有单元端点n1个：把边界2离散为n2个 单元，有单元端点n2个，则边界「共划分为n (n=m,+n2)个单元，共有n个单元端点。在实践 中，边界单元一般被离散为常单元、线性单元或者 高次单元，也有采用样条插值方式的边界单元。 在本文中，采用线性单元就二维问题进行计算。 图1研究区城 基函数选取如下： 设想包含在Ω内的不可压缩粘性流体整体嵌 [L+E,x∈T-i 2 入在平面无限域上的不可压缩粘性流体之中。本 P= 文只计算有界域上的流动，故假定在无穷远处流 1- ,xeT,其中ee[-1,1],i=1,…,n 2 速为零。由参考文献的，经计算可得到基于单层 0,其他 位势的边界积分方程： (y)= 名.x-(pden2 设1=含4=含”e其中：0为 t,在节点i处的值。即得离散形式： 其中t=(t1,t2)是待定的密度函数，其分量t:是 含0，(u,p穿超「时的跃变，类似地可推出流 体压力场的积分表达式，但本文约去压力场的 ds,+ uizit ds,)=b + 计算。 2变分方程的形成及处理 Jr-i+「r 2.1形成变分方程 (5) 由单层位势出发来研究边界量u。和边界量t 其中： ?1994-2017 China Academic Journal Electronic Publishing House.All rights reserved.http://www.cnki.net
意点的流速值。本文主要讨论其中用利用 Fortran 计算所求点的流速与用 Matlab 画出相应的流线图 的主要过程。 1 复连通区域 Stokes 问题的边界积分方程 设 Γ1 是 R2 中有限区域 Ω 的封闭曲线，Γ2 是 Γ1 形成的封闭区域中的一条封闭曲线，Γ2 形成的 区域为 Ω + 0 ，Ω － 0 = Ω － Ω + 0 。设 Γ = Γ1 + Γ2 ，考虑如 下的问题: － υΔu + gradp = 0， 在 Ω 中 div u = 0， 在 Ω 中 u Γ1 +Γ2 = u0， 在 Γ1 + Γ { 2 上 ( 1) 其中: u = ( ) u1，u2 是流速; p 是压强; u0 是闭曲线边 界 Γ1 以及内部闭曲线边界 Γ2 上的已知函数。 图 1 研究区域 设想包含在 Ω 内的不可压缩粘性流体整体嵌 入在平面无限域上的不可压缩粘性流体之中。本 文只计算有界域上的流动，故假定在无穷远处流 速为零。由参考文献［5］，经计算可得到基于单层 位势的边界积分方程: uk ( ) y = ∑ 2 i = 1 ∫Γ Uik ( ) x － y ti ( ) x dsx，y ∈ Ω ( 2) 其中 t = t ( ) 1，t2 是待定的密度函数，其分量 ti 是 ∑ 2 j = 1 σij( ) u，p nj 穿越 Γ 时的跃变。类似地可推出流 体压力场的积分表达式，但本文约去压力场的 计算。 2 变分方程的形成及处理 2． 1 形成变分方程 由单层位势出发来研究边界量 u0 和边界量 t 的关系，由于 u 在穿越边界时的连续性，从式( 2) 得到联 系 已 知 函 数 u0 和 中 间 变 量 t 的 积 分 方 程组: u0k ( y) = ∑ 2 i = 1 ∫Γ Uik ( x － y) ti ( x) dsy， k = 1，2，y ∈ Γ ( 3) 对于方程( 3) ，方程两边同乘以 t'( ) y ，然后在 Γ 上 积分，便得到变分方程组: ∑ 2 i，k = 1 ∫Γ t' k ( ) y ∫Γ Uki ( x，y) ti ( x) dsxdsy = ∑ 2 k = 1 ∫Γ U0k ( y) t' k ( y) dsy ( 4) 2． 2 离散变分方程 为了数值求解变分方程( 4) ，边界 Γ 必须离 散为一系列单元，从而把边界积分方程转化为线 性代数方程组进行求解。设把边界 Γ1 离散为 n1 个单元，有单元端点 n1 个; 把边界 Γ2 离散为 n2 个 单元，有单元端点 n2 个，则边界 Γ 共划分 为 n ( ) n = n1 + n2 个单元，共有 n 个单元端点。在实践 中，边界单元一般被离散为常单元、线性单元或者 高次单元，也有采用样条插值方式的边界单元。 在本文中，采用线性单元就二维问题进行计算。 基函数选取如下: φi = 1 + ε 2 ，x ∈ Γi－1 1 － ε 2 ，x ∈ Γi，其中 ε ∈［－ 1，1］，i = 1，…，n 0， 其他 设 t1 = ∑ n i = 1 t ( i) 1 φi，t2 = ∑ n i = 1 t ( i) 2 φi，其中: t ( i) 1 、t ( i) 2 为 t1，t2 在节点 i 处的值。即得离散形式: ∑ n i = 1 ∫Γj－1 +Γj φj ( ) y dsy * ∫Γi－1 +Γi u* 11φi t ( i) 1 dsx + ∫Γi－1 +Γi u* 12φi t ( i) ( ) 2 dsx = b( j) 1 ∑ n i = 1 ∫Γj－1 +Γj φj ( ) y dsy * ∫Γi－1 +Γi u* 21φi t ( i) 1 dsx + ∫Γi－1 +Γi u* 22φi t ( i) ( ) 2 dsx = b( j) 2 ( 5) 其中: 118 重庆理工大学学报
王小军，等：复连通二维Stokes问题的Galerkin边界元解法 119 的坐标以及所求解。 在输入功能子程序NPUT中要实现的功能是 输入边界节点坐标、已知边界条件和所需计算的 j=1,2,…,n1,n1+1,…,n 内点坐标。 把上述线性方程组写成简洁形式：AU=B,其中 在形成系数矩阵的子程序FMAT中，首先通 过两层循环采用4点Gauss数值积分公式形成右 A =[ayl nxem'U= 端项，然后通过3层循环调用解析积分公式和4 点Gauss数值积分公式形成系数矩阵，其中在第1 =1,2,…,n 层积分中，存在奇异性时调用解析积分公式，不存 在奇异性时调用4点Gauss数值积分公式。 2.3变分方程中重积分的处理 计算解析积分公式的子程序FUN中，通过判 要求解上述线性方程组(5)，最主要的问题是 断语句根据积分源点到有向线段的距离是否为零 求解系数矩阵。然而由于积分存在奇异性，若采 分情况计算系数矩阵对应元素的数值。 用数值积分，必然产生较大的误差。为了避免出 计算4点Gauss数值积分公式的子程序，调用 现这样的问题，在求二重积分时，第1重积分存在 为Gauss数值积分作准备的子程序，将坐标转变为 奇异性时采用精确解析积分，不存在奇异性时以 无因次积分点坐标，实现非奇异情况下的数值 及第2重积分采用Gauss数值积分。在利用Gauss 积分。 数值积分计算第2重积分时，对于以闭曲线为边 计算线性代数方程组的源程序，采用列主元 界的区域问题而言，基函数无奇异性，可以方便的 消去法求解边界积分方程组。 使用通常的Gauss数值积分公式。解析公式的推 计算流速的子程序中采用2层循环，将求得 导方法与具体公式见参考文献2]。利用边界积 的密度数值带入解的表达式求得未知点的流速及 分方程得到密度函数后，将其数值带入解的表达 其方向。 式即得流速。 输出结果的子程序将计算出的结果写入文本 文件，以便对比和画流线图。 3程序设计 Matlab程序中首先对区域进行剖分，得到所 求点坐标，将坐标输出到文本文件，由Fortran程序 鉴于不同程序设计语言的优势，程序采用两 计算出流速后，加载并读取保存结果的文本文件， 种语言，首先用Fortran计算出所求点的流速，保存 然后调用已知函数画出流线图，最后画出边界并 结果在文本文件中，然后用Matlab读取画出流 确定坐标名称。 线图。 Fortran程序采用模块化设计，包括主程序M- 4数值算例 AN以及8个子程序，其中主程序控制程序的执行 顺序，采用多个输入输出文本文件读取写入以便 算例1单位圆流速为(0,0)，矩形区域 于与Matlab结合。在主程序中设置全局变量，包 [-5,5,-10,10]边界上的速度为(1,0)，计算矩 括了边界点，边界点对应的流速，边界单元长度， 形内单位圆以外的点的流速并画出流线图。 稀疏矩阵，右端项，单层位势的向量密度，所求点 在此算例中，内部圆边界划分为16个单元， ?1994-2017 China Academic Journal Electronic Publishing House.All rights reserved.http://www.cnki.net
bj = b( j) 1 b( ) ( j) 2 = ∫Γj－1 +Γj u01φj dsy ∫Γj－1 +Γj u02φj ds y ， j = 1，2，…，n1，n1 + 1，…，n 把上述线性方程组写成简洁形式: AU = B，其中: A = [ ] aij ( 2n) ×( 2n) ，U = t ( i) 1 t [ ] ( i) 2 b = b( i) 1 b [ ] ( i) 2 ，i = 1，2，…，n 2． 3 变分方程中重积分的处理 要求解上述线性方程组( 5) ，最主要的问题是 求解系数矩阵。然而由于积分存在奇异性，若采 用数值积分，必然产生较大的误差。为了避免出 现这样的问题，在求二重积分时，第 1 重积分存在 奇异性时采用精确解析积分，不存在奇异性时以 及第 2 重积分采用 Gauss 数值积分。在利用 Gauss 数值积分计算第 2 重积分时，对于以闭曲线为边 界的区域问题而言，基函数无奇异性，可以方便的 使用通常的 Gauss 数值积分公式。解析公式的推 导方法与具体公式见参考文献［2］。利用边界积 分方程得到密度函数后，将其数值带入解的表达 式即得流速。 3 程序设计 鉴于不同程序设计语言的优势，程序采用两 种语言，首先用 Fortran 计算出所求点的流速，保存 结果在文本文件中，然后用 Matlab 读取画出流 线图。 Fortran 程序采用模块化设计，包括主程序 MIAN 以及 8 个子程序，其中主程序控制程序的执行 顺序，采用多个输入输出文本文件读取写入以便 于与 Matlab 结合。在主程序中设置全局变量，包 括了边界点，边界点对应的流速，边界单元长度， 稀疏矩阵，右端项，单层位势的向量密度，所求点 的坐标以及所求解。 在输入功能子程序 INPUT 中要实现的功能是 输入边界节点坐标、已知边界条件和所需计算的 内点坐标。 在形成系数矩阵的子程序 FMAT 中，首先通 过两层循环采用 4 点 Gauss 数值积分公式形成右 端项，然后通过 3 层循环调用解析积分公式和 4 点 Gauss 数值积分公式形成系数矩阵，其中在第 1 层积分中，存在奇异性时调用解析积分公式，不存 在奇异性时调用 4 点 Gauss 数值积分公式。 计算解析积分公式的子程序 FUN 中，通过判 断语句根据积分源点到有向线段的距离是否为零 分情况计算系数矩阵对应元素的数值。 计算 4 点 Gauss 数值积分公式的子程序，调用 为 Gauss 数值积分作准备的子程序，将坐标转变为 无因次 积 分 点 坐 标，实现非奇异情况下的数值 积分。 计算线性代数方程组的源程序，采用列主元 消去法求解边界积分方程组。 计算流速的子程序中采用 2 层循环，将求得 的密度数值带入解的表达式求得未知点的流速及 其方向。 输出结果的子程序将计算出的结果写入文本 文件，以便对比和画流线图。 Matlab 程序中首先对区域进行剖分，得到所 求点坐标，将坐标输出到文本文件，由 Fortran 程序 计算出流速后，加载并读取保存结果的文本文件， 然后调用已知函数画出流线图，最后画出边界并 确定坐标名称。 4 数值算例 算例 1 单 位 圆 流 速 为 ( 0，0 ) ，矩 形 区 域 ［－ 5，5，－ 10，10］边界上的速度为( ) 1，0 ，计算矩 形内单位圆以外的点的流速并画出流线图。 在此算例中，内部圆边界划分为 16 个单元， 王小军，等: 复连通二维 Stokes 问题的 Galerkin 边界元解法 119
120 重庆理工大学学报 外部闭曲线划分为60个单元，边界总剖分为76 r(1,0) -10≤x≤10，y=5 (u,)= 单元，由Matlab剖分在矩形区域内确定1521个 (0,0)x=±10，ory=-5 点，除去圆内部余1492个有效点，由Fortran程序 计算矩形内单位圆以外的点的流速并画出流 计算出流速后输出到文本文件，耗费时间以秒为 线图。 单位，再由Matlab程序读取画出流线图，耗费时间 在此算例中，边界剖分与计算流程与上述两 也是以秒为单位，与其他方法计算结果相对比精 例完全一致，计算、画图耗用时间以及结果精确度 度较高，画出的流线图与实际情况符合很 也完全一致。不同点在于，Matlab画图中放大后 好（图2）。 将清晰的出现2个小的漩涡。漩涡如图4所示。 42924.10 图2算例1流线图 算例2单位圆流速为(0,0)，矩形区域 [-5,5,-10,10]边界上的速度为 r(0,0) y=±5 25 (u,)= 20 (1-y2/25,0)x=±10 15 计算矩形内单位圆以外的点的流速并画出流 10 线图。 在此算例中，边界剖分与计算流程与算例1 完全一致，计算、画图耗用时间以及结果精确度也 完全一致（图3）。 15 图4算例3的漩涡图 0 2468 通过以上数值算例可以看出，该种方法精度 图3算例2流线图 较高，计算结果符合客观事实。 算例3单位圆流速为(0,0)，矩形区域 [-5,5,-10,10]边界上的速度为 (下转第126页) ?1994-2017 China Academic Journal Electronic Publishing House.All rights reserved.http://www.cnki.net
外部闭曲线划分为 60 个单元，边界总剖分为 76 单元，由 Matlab 剖分在矩形区域内确定 1 521个 点，除去圆内部余 1 492个有效点，由 Fortran 程序 计算出流速后输出到文本文件，耗费时间以秒为 单位，再由 Matlab 程序读取画出流线图，耗费时间 也是以秒为单位，与其他方法计算结果相对比精 度较 高，画出的流线图与实际情况符合很 好( 图 2) 。 图 2 算例 1 流线图 算 例 2 单 位 圆 流 速 为 ( 0，0 ) ，矩 形 区 域 ［－ 5，5，－ 10，10］边界上的速度为 ( ) u，v = ( ) 0，0 y = ± 5 1 － y 2 ( ) /25，0 x = ± { 10 计算矩 形 内 单 位 圆 以 外 的 点 的 流 速 并 画 出 流 线图。 在此算例中，边界剖分与计算流程与算例 1 完全一致，计算、画图耗用时间以及结果精确度也 完全一致( 图 3) 。 图 3 算例 2 流线图 算 例 3 单 位 圆 流 速 为 ( 0，0 ) ，矩 形 区 域 ［－ 5，5，－ 10，10］边界上的速度为 ( ) u，v = ( ) 1，0 － 10 ≤ x ≤ 10，y = 5 ( ) 0，0 x = ± 10，or y = － { 5 计算矩 形 内 单 位 圆 以 外 的 点 的 流 速 并 画 出 流 线图。 在此算例中，边界剖分与计算流程与上述两 例完全一致，计算、画图耗用时间以及结果精确度 也完全一致。不同点在于，Matlab 画图中放大后 将清晰的出现 2 个小的漩涡。漩涡如图 4 所示。 图 4 算例 3 的漩涡图 通过以上数值算例可以看出，该种方法精度 较高，计算结果符合客观事实。 ( 下转第 126 页) 120 重庆理工大学学报
126 重庆理工大学学报 振动系统的振幅与初始位置和初始速度有关系。初始位置增大，振幅增大；初始位置不变，若初始速度增 大，则振幅也增大。 2)在其x方向和在y方向的振动曲线都可以看成是一个形变了的余弦曲线，其上附加有依赖于振 幅的高频颤动。 3结束语 对弹簧系统中常见的二维非线性振动问题，可利用拉格朗日方法得到其振动控制微分方程，借助于 计算机和maple语言在计算方面的超强功能，成功解决了该类非线性振动问题。这种方法有效简便，这 就为非线性问题探索出了一种极好的求解途径。 振动曲线彼此相似，其波形与振幅无关，具体形状介于余弦波和三角波之间，可以看成是一个变形了 的余弦波，其上附加有依赖于振幅的高频颤动。 参考文献： 0]廖旭，任学藻.组合线性弹簧振子中的非线性振动叮.大学物理，2008,27(2)：25-28. D]包兴明，周志坚，袁玉全.弹簧系统一种常见的二维非线性振动们.西南师范大学学报，2008,33(4)：28-31. B]卓崇培.非线性物理学0.天津：天津科学技术出版社，1996：26-61. 4]刘式适，刘式达.物理学中的非线性方程0.北京：北京大学出版社，2000：6,66 [5]黎捷.MAPLE9.0符号处理及应用M.北京：科学出版社，2004. (责任编辑刘射) (上接第120页) 类边界积分方程0].计算数学，2005,2(1)：1-10. )祝家麟.椭圆边值问题的边界元分析门.北京：科 参考文献： 学出版社，1991. 5]祝家麟.定常toks问题的边界积分方程法们.计 [Hsiao G C.A modified Galerkin scheme for equations 算数学，1986(8)：281-289. with natural boundary conditions [C]//Vichnevetsky R, 6]Zhu Jialin.The boundary integral equation method for Vigness J.Numerical Mathematics and Applications.B. solving stationary Stokes problems [C]//Feng Kang. V:Elsevier Science Publishers,1986:193-197. Proc.of the 1984 Beijing Symp.On Diff.Geometry and [D]向瑞银.平面定常Stokes方程的Galerkin边界元解 Diff.Equations.Beijing:Science Press,1985. 法].重庆大学学报，2006(2)：29-30. B]张耀明，温卫东.平面定常Stokes问题的无奇异第一 (责任编辑刘舸) ?1994-2017 China Academic Journal Electronic Publishing House.All rights reserved.http://www.cnki.net
振动系统的振幅与初始位置和初始速度有关系。初始位置增大，振幅增大; 初始位置不变，若初始速度增 大，则振幅也增大。 2) 在其 x 方向和在 y 方向的振动曲线都可以看成是一个形变了的余弦曲线，其上附加有依赖于振 幅的高频颤动。 3 结束语 对弹簧系统中常见的二维非线性振动问题，可利用拉格朗日方法得到其振动控制微分方程，借助于 计算机和 maple 语言在计算方面的超强功能，成功解决了该类非线性振动问题。这种方法有效简便，这 就为非线性问题探索出了一种极好的求解途径。 振动曲线彼此相似，其波形与振幅无关，具体形状介于余弦波和三角波之间，可以看成是一个变形了 的余弦波，其上附加有依赖于振幅的高频颤动。 参考文献: ［1］ 廖旭，任学藻． 组合线性弹簧振子中的非线性振动［J］． 大学物理，2008，27( 2) : 25 － 28． ［2］ 包兴明，周志坚，袁玉全． 弹簧系统一种常见的二维非线性振动［J］． 西南师范大学学报，2008，33( 4) : 28 － 31． ［3］ 卓崇培． 非线性物理学［M］． 天津: 天津科学技术出版社，1996: 26 － 61． ［4］ 刘式适，刘式达． 物理学中的非线性方程［M］． 北京: 北京大学出版社，2000: 6，66． ［5］ 黎捷． MAPLE 9． 0 符号处理及应用［M］． 北京: 科学出版社，2004． ( 责任编辑 刘 舸 檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶檶 ) ( 上接第 120 页) 参考文献: ［1］ Hsiao G C． A modified Galerkin scheme for equations with natural boundary conditions［C］/ /Vichnevetsky R， Vigness J． Numerical Mathematics and Applications． B． V: Elsevier Science Publishers，1986: 193 － 197． ［2］ 向瑞银． 平面定常 Stokes 方程的 Galerkin 边界元解 法［J］． 重庆大学学报，2006( 2) : 29 － 30． ［3］ 张耀明，温卫东． 平面定常 Stokes 问题的无奇异第一 类边界积分方程［J ］． 计算数学，2005，2( 1) : 1 － 10． ［4］ 祝家麟． 椭圆边值问题的边界元分析［M］． 北京: 科 学出版社，1991． ［5］ 祝家麟． 定常 Stokes 问题的边界积分方程法［J］． 计 算数学，1986( 8) : 281 － 289． ［6］ Zhu Jialin． The boundary integral equation method for solving stationary Stokes problems［C］ / /Feng Kang． Proc． of the 1984 Beijing Symp． On Diff． Geometry and Diff． Equations． Beijing: Science Press，1985． ( 责任编辑 刘 舸) 126 重庆理工大学学报