轴向均布荷载压杆稳定问题的重心插值配点法

重心插值配点法及其应用摘要：重心Lagrange插值具有数值稳定性好、计算精度高的优点。

采用重心Lagrange插值近似未知函数，建立未知函数各阶导数的微分矩阵。

采用微分矩阵近似未知函数的导数，利用配点法将控制方程和边界条件离散为代数方程组，通过求解代数方程组，从而可求解偏微分方程。

数值算例表明，重心插值配点法具有原理简单，易于程序实现和数值计算精度高的优点。

关键词：重心Lagrange插值；微分矩阵；配点法Barycentric Lagrange interpolation collocation Method and its ApplicationAbstract：Barycentric Lagrange interpolation collocation method has excellent numerical stability and high accuracy. this paper presents the Barycentric Lagrange interpolation collocation method to get the differentiation Matrix of unknown function. So the control equation can be expressed as linear systems by the collocation method. According to those formulas, differential equations can be soluted. The principle of this method is simple and easy to programming. The accuracy and the numerical stability are very excellent.Key words: Barycentric Lagrange interpolation, differentiation Matrix, collocation method,0 引言具有初、边值条件的常、偏微分方程的解析解常无法通过理论推导获得，一种有效的途径是采用数值求解方法获得具有一定数值精度的近似解，这些数值方法包括：有限差分法、有限单元法、边界单元法、无网格法及一些杂交使用的方法等数值求解方法。

梁方程降阶计算的重心插值配点法徐子康;王兆清;孙浩森;李金【摘要】采用重心插值配点法求解梁方程时,随着计算节点数量的持续增加,其计算精度将逐步下降.通过对降阶计算重心插值配点法的研究,可为数值求解梁方程提供一种数值稳定性好、计算精度高的新方法.文章基于重心Lagrange插值及其微分矩阵,推导了梁方程降阶计算重心插值配点法的公式,并通过数值算例验证其有效性.结果表明:随着计算节点数量的持续增加,降阶法的计算精度仍保持在10-10 ～ 10-12范围内;求解两端简支的梁方程时,两步降阶法的计算精度高于一步降阶法;直接法计算矩阵条件数与节点数的7次方是同阶的,而一步降阶法计算矩阵条件数与节点数的4次方是同阶的,降阶法可以有效地降低计算矩阵的条件数,提高计算精度;重心插值配点法采用矩阵—向量形式的计算公式,便于程序的编写,提高了计算效率.%When the beam equations are solved by barycentric interpolation collocation method,computational accuracy will decline gradually as the number of nodes increases.The research on barycentric interpolation collocation method based on depression of order,can provide the new method that has a good numerical stability and high computational accuracy for beam equations.Based on barycentric Lagrange interpolation and its differential matrices,the formula of barycentric interpolation collocation method based on depression of order is derived.Numerical examples are given to verify the effectiveness of the proposed method.The result shows that the computational accuracy of depression of order method remains the range of 10-10 ～ 10-12 as the number of nodes increases.When the beam equations with simply supported ends aresolved,the computational accuracy of the two-step depression of order method is higher than the one-step depression of order method.The condition number of the direct method is close to 7th power of the number of nodes,and the condition number of one-step depression of order method is close to 4th power of the number of nodes.The depression of order method can effectively reduce the condition number of computing matrix such that improves the computational accuracy.By applying the computational formula with matrix-vector form,the program is easy to write and the computational efficiency of barycentric interpolation collocation method can be improved remarkably.【期刊名称】《山东建筑大学学报》【年(卷),期】2017(032)003【总页数】6页(P245-250)【关键词】梁方程;降阶法;重心Lagrange插值;配点法【作者】徐子康;王兆清;孙浩森;李金【作者单位】山东建筑大学工程力学研究所,山东济南250101;山东建筑大学工程力学研究所,山东济南250101;山东建筑大学学报编辑部,山东济南250101;山东建筑大学理学院,山东济南250101【正文语种】中文【中图分类】O302梁方程是工程技术领域中常见的四阶常微分方程。

含端部弹性约束和铰支约束压杆的稳定性问题研究雷明伟;骆凯;史文谱【摘要】弹性压杆在工程中有广泛应用,压杆的失稳问题是压杆失效的重要原因之一,是压杆可靠性优化设计中需要重点关注的问题,临界载荷的计算方法和表达式是数学建模和优化的基础.本文针对含端部弹性约束和铰链支撑约束压杆的临界失稳问题进行探讨,提出了解析分析方法,得到的结果当弹性约束刚度系数分别为无穷大和零时,将分别简化为一端固定一端铰支和两端铰支的压杆临界失稳载荷计算结果,本文的分析方法和结论对含多个弹性约束的压杆失稳问题研究都有一定的参考意义.【期刊名称】《烟台大学学报（自然科学与工程版）》【年(卷),期】2018(031)002【总页数】5页(P153-157)【关键词】压杆稳定;弹性约束;临界载荷【作者】雷明伟;骆凯;史文谱【作者单位】烟台大学机电汽车工程学院,山东烟台264005;烟台大学机电汽车工程学院,山东烟台264005;烟台大学机电汽车工程学院,山东烟台264005【正文语种】中文【中图分类】TB12压杆稳定问题是机械工程中广泛存在的问题,许多工程失效事故常常是由于压杆抗失稳能力不够造成的,比如大型桁架结构的局部坍塌、桥梁断裂事故、油压千斤顶、螺杆传动的压榨机构、含有连杆的传动机构等．文献[1]针对Q500qE高强度钢压杆稳定问题进行探讨,利用有限元软件Abaqus建立实体和壳单元压杆计算模型,数值分析中考虑了非线性、初始缺陷、焊接应力、杆件长细比及杆件截面形状等因素;文献[2]利用有限元软件ANSYS分析了初始挠度及中间弹性支承对压杆临界失稳载荷的影响问题;文献[3]针对常规有限元法不能得到稳定问题精确解的情况提出了新的压杆稳定单元,建立了精确单元形函数,利用迭代算法确定临界载荷和失稳模态．事实上,确定压杆临界载荷的方法除了上述有限元等数值计算方法外,还有微分方程求解法[4]、积分方程法[5]、微分变换法[6]、泰勒级数法[7]、重心插值配点法[8]、刚度判别法[9]以及可为压杆选用型材的图解算法[10]等．从上述列出文献以及其他未列出的参考文献来看,等截面压杆稳定性问题解析求解的难点在于压杆端部约束的类型和性质,目前两端固定、两端铰支、一端固定一端铰支、一端固定一端自由的压杆临界失稳载荷已经给出解析解[4],但有关含有端部弹性约束的压杆稳定性问题的解析研究却未见有论文发表,而实际上端部固定约束只是一种理想的力学模型,这是一种刚性约束,现实中是不存在的．为此本文采用近似挠度微分方程,通过设定合理的约束边界条件,利用求解微分方程的办法得到了确定压杆临界载荷的超越方程,再利用迭代的办法得到问题的解．当假设弹性约束刚度系数很大或很小时,问题将分别简化为人们熟知的一端固定一端铰支的压杆稳定性问题和两端铰支的压杆稳定性问题,得出的结论与已知结果是严格一致的,说明本文分析方法和结论是正确可行的,本文方法和研究思路对于其他类似约束性质的压杆稳定问题研究都具有理论参考意义．1 问题模型及理论分析如图1所示,一根均质等径压杆,长度为L,材料杨氏模量为E,截面惯性矩为I．其一端O处为弹性支撑(弹性刚度系数为K),另外一端A处为铰链支撑,受到图示轴向压力F的作用,发生如图1所示的弯曲变形．建立图示坐标系xOy,其中x轴沿着压杆原始位置时的轴心线方向,y轴垂直于压杆轴心线方向．按照材料力学分析方法以及有关物理量的符号规定,压杆任意位置x处的截面弯矩M(x)可表为M(x)=-Fw+T(x-L),(1)其中:w为压杆失稳弯曲变形挠度,T是由于压杆受到压缩作用失稳发生弯曲导致左端弹性约束处产生的抵抗弯矩M0的作用下为保证压杆静力平衡的右端支撑横向作用力．并假设F>0,T>0,即取两者的绝对值．图1 一端弹性约束一端铰支压杆的失稳变形Fig.1 Unstable deformation of buckling column with elasticside and pin support压杆挠度满足的近似微分方程为EIw″=M(x).(2)将方程(1)代入(2)中整理有EIw″+Fw=T(x-L),(3)或者w″+ FEI w= FEI (x-L).根据常微分方程理论,其解为通解w1(x)和特解w2(x)之和．通解w1(x)即为方程(3)对应的齐次常微分方程的解,可表为w1(x)=c1coskx+c2sinkx,(4)其中:k= FEI ,cj(j=1,2)为待定系数．按照非齐次常微分方程的求解方法和理论,根据现有问题的特点,可假设其特解为w2(x)=ax+b,(5)其中:a,b是2个待定系数．将方程(5)带入微分方程(3)中可得待定系数的解为a=T/F ,b=-TL/F .(6)这样微分方程(3)的全解为w(x)=c1coskx+c2sinkx+Tx/F-TL/F.(7)从图1可看出,问题的边界条件可表为w(0)=0 ,w(0)=-α ,w(L)=0 ,(8)其中:α是压杆左端面外法线转过的角度大小．此外根据压杆的静力平衡条件有Kα=T/L,(9)由此得T=Kα/L.(10)将方程(7)代入边界条件(8)中的第一个方程有c1=TL/F .此外,将方程(7)两边微分有w′(x)=-kc1sinkx+kc2coskx+T/F.(11)将方程(11)代入边界条件(8)中第2个边界条件有c2=(-T/F-α)/k.将方程(7)代入边界条件(8)中的第3个边界条件可得tg(kL)=-c1/c2= kL1+FL/K .(12)显然,作为特例,也是验证本文推导结论是否正确的一种方法,当K=∞时,即压杆左端为固定,这样原有的问题成为一端固定一端铰支的压杆问题,而该问题的结论是有的,结果是tg(kL)=kL.(13)当K=0时,即压杆左端也是铰链支撑,则有tg(kL)=0或者sin(kL)=0.(14)从方程(13)和方程(14)可看出,本文分析的问题的简化结果与文献[4]所给结论是一致的．为了下面讨论方便起见,假设β=kL,重新整理方程(12)有tg β= β1+EIβ2/(LK) =f(β,K).(15)2 算例结果及分析作为数值算例,假设压杆为45号钢制成,其比例极限为σp=280 MPa;E=210GPa;L=1.5 m,d=30 mm.欧拉公式可以应用的柔度下限为λ1=π E/σp =86.0,压杆约束介于一端固定一端铰支和两端铰支的情形之间,故有效长度系数(或长度因数)为μ∈(0.7,1.0),压杆截面为圆形,其截面惯性半径:L= I/S =d/4=7.5 mm ,其中:S是压杆横截面积．截面惯性矩:I=i2S=πd4/64=3.974×10-8 m4,EI=210×109×3.974×10-8=8 345.4 N·m2.计算柔度:λmin=μL/i=0.7×1.5/0.007 5=140>λ1.分别取:K=1,10,20 N·m·rad-1,K=100,200,300 N·m·rad-1,K=1 000,2 000,3 000 N·m·rad-1,K=104,2×104,3×104N·m·rad-1,K=105,2×105,3×105N·m·rad-1,K=106,2×106,3×106N·m·rad-1,计算结果分别如图2～7所示．在这些计算结果中,切线斜率最大的曲线是函数y=tg β的图像,其他曲线是函数y=f(β,K)的图像．从结果来看,随着K的增大,函数曲线从下向上依次排开,由于每一组K值仅仅取了3个数值,故每个图中都有3条曲线(不算正切曲线)．由于正切函数是周期函数,其最小正周期为π,在正实数轴范围内其间断点为jπ/2(j=1,3,5,7,…)．它在每个周期范围内均为单调增函数．从图2来看,由于刚度系数K较小,正切函数曲线和这些不同K值的函数f(β,K)曲线的交点都非常接近于零,这是一个平凡解,没有实际意义;但是根据压杆失稳的物理特性,压杆临界失稳力是客观存在的,所以只能从sinβ=0方程中得到解答;为了进一步说明问题,又给出了计算结果图3和图4,在β∈[0,π/2]范围内,情况与图1是类似的,只是函数f(β,K)的大小有所不一样,显然随着K的增大,函数f(β,K)的曲线呈现向上移动的趋势,但能够看出在K≤3 000 N·m·rad-1范围内,正切曲线y=tg β与函数f(β,K)在β∈[0,π]范围内除了坐标原点附近外不可能出现其他交点了,因此只能借助于方程sinβ=0来确定压杆的临界力了,即对于本文讨论的算例来说,K≤3 000 N·m/rad范围内,弹性约束端可看作铰链支撑约束了,即原来的一端弹性约束一端铰链约束的压杆稳定性问题可按照两端铰支约束压杆问题来处理．图2 函数y=tg β与函数y=f(β,K)的相交Fig.2 Intersection of function y=tg β and y=f(β,K) with dif-ferent图3 函数y=tg β与函数y=f(β,K)的相交Fig.3 Intersection of function y=tg β and y=f(β,K) with dif-ferent图4 函数y=tg β与函数y=f(β,K)的相交Fig.4 Intersection of function y=tg β and y=f(β,K) with different当弹性刚度系数K≥104 N·m·rad-1时,结果如图5～7所示,函数f(β,K)与正切函数tg β在[π,3π/2)范围内有非零交点．容易看出,随着刚度系数K的增大,这个交点的横坐标和纵坐标都随着增大,即该交点向右上方移动,并且当K越来越大时,正切函数tg β曲线与函数f(β,K)的曲线的交点逐步稳定下来,其横坐标越来越接近于4.49,这说明当K增大到一定程度时,压杆原来的弹性约束端逐步演变为固定端约束了,这与前面的理论分析是一致的,而且容易看出来,当刚度系数K增大到一定程度时,函数f(β,K)就趋近于斜率为1的直线了,这与事实和理论分析是吻合的(如图7所示)．图5 函数y=tgβ与函数y=f(β,K)的相交Fig.5 Intersection of function y=tgβ and y=f(β,K) with dif-ferent图6 函数y=tgβ与函数y=f(β,K) 的相交Fig.6 Intersection of function y=tgβ and y=f(β,K) with dif-ferent图7 函数y=tg β与函数y=f(β,K) 的相交Fig.7 Intersection of function y=tg β and y=f(β,K) with dif-ferent3 结论从前面的理论分析和数值算例来看,可得如下几个结论:(1)在弹性刚度系数K较大的范围内(比如对于本文算例来说,K≤3 000 N·m·rad),压杆弹性约束端可近似简化为铰链支撑约束,因而可利用两端铰支压杆的临界力公式[4]处理压杆失稳问题;(2)弹性约束刚度系数K对失稳临界力是有影响的,只是不特别敏感而已;随着刚度系数K的增大,临界失稳力是单调增加的;(3)当刚度系数K增大到一定程度时,它对压杆失稳临界力的影响逐步减弱,直到稳定下来(比如对于本文讨论的算例而言,K≥3×106 N·m·rad),K(β,K)的影响可以忽略了,使得K(β,K)→β．即原来的一端弹性约束一端铰链支撑约束压杆失稳问题可以高精度地近似看作一端固定一端铰链支撑约束压杆的稳定性问题了．(4)刚度系数K对于本文讨论的压杆失稳问题的临界载荷客观上还是有影响的,当这种影响不可忽略时,完全可按照本文得到的结果,进一步利用数值迭代方法得到相应失稳临界力的解答．参考文献:[1] 鞠晓臣, 田越, 赵欣,等. Q500qE高强钢压杆稳定研究[J]. 铁道建筑,2015(10):80-84.[2] 张晓霞, 钟文生, 姚远. 初始挠度及中间弹性支承对压杆稳定的影响分析[J]. 设计与研究, 2011,38(6): 1-4.[3] 任风鸣, 范学明. 弹性压杆稳定问题的精确解法[J]. 建筑科学, 2008, 24(3): 12-14.[4] 刘鸿文. 材料力学[M].第5版. 北京: 高等教育出版社, 2011, 290-303.[5] 陈春. 积分方程在压杆稳定中的应用[J]. 重庆建筑工程学院学报, 1991, 13(4): 74-78.[6] 禹金云, 胡辉. 微分变换在压杆稳定问题中的应用[J]. 湘潭师范学院学报(自然科学版), 2004,26(1): 38-40.[7] 张适, 谢冬梅. 求压杆稳定的临界力的数学方法[J]. 云南民族学院学报(自然科学版), 1996, 5(1): 25-30.[8] 于卫涛, 宋洁, 赵维霞. 轴向均布荷载压杆稳定问题的重心插值配点法[J]. 山东建筑大学学报, 2011,26(4): 353-355.[9] 周奇才, 李文军, 吴青龙,等. 空间压杆极值失稳轴向刚度判别法[J]. 机械强度,2016, 38(1):94-98.[10] 焦良. 为压杆选用型材的图解法[J]. 机械强度,1992,14(1): 64-65.。

考虑轴向均布载荷时压杆的稳定性计算黄开志;陈小亮【摘要】为了求得压杆同时承受轴向均布载荷和集中载荷时,临界载荷的计算公式,首先对仅承受轴向均布载荷的压杆,用初参数法,导出了临界载荷特征方程,由软件分析特征方程发现,"固支-定向"、"固支-自由"、"铰支-定向"支承的压杆,轴向均布载荷对其稳定性有明显影响,并求得了临界载荷的近似解;其次采用载荷换算与叠加的方法,求得了压杆同时承受轴向均布载荷和集中载荷时,临界载荷计算的经验公式;最后就"固支-自由"支承的压杆,与其他一些研究结果进行了比较,本文结果与"平均结果"较吻合.【期刊名称】《力学与实践》【年(卷),期】2015(037)006【总页数】6页(P744-749)【关键词】材料力学;压杆;稳定性;临界载荷;经验公式;欧拉公式【作者】黄开志;陈小亮【作者单位】重庆科技学院数理学院,重庆 401331;重庆科技学院数理学院,重庆401331【正文语种】中文【中图分类】O341对同时承受轴向均布载荷和集中载荷的压杆稳定性问题，文献［1］指出，难以求得临界载荷的精确解，该文作者采用摄动法得到了临界载荷的特征方程，给出了临界载荷的一些数值解；文献［2］首先假设了压杆挠曲线近似为正弦函数，然后用能量原理得到了临界载荷的解答；文献［3］用幂级数法求得了临界载荷的解答.在这些近似方法中，较详细研究了“固支--自由”支承的情况，对其他支承情况，没有提供详实的临界载荷公式.本文拟先导出压杆仅承受轴向均布载荷q时，其临界载荷的近似解，并根据不同的支承情况作载荷换算，再与压杆轴向集中载荷叠加，以求得压杆同时承受轴向均布载荷和集中载荷时，其临界载荷计算的经验公式.设长为l、抗弯刚度为E I的压杆，仅承受轴向均布载荷q时处于微弯曲平衡状态，其受力和变形情况可简化为图1所示模型.1.1 变形方程在微弯曲平衡状态的压杆中取d x微段，因小变形，可视d x段为直线，按材料力学习惯设内力为正，其受力情况如图2所示，则静力平衡方程为因小变形，故sin w′（x）≈w′（x），并略去二阶微量，上述二式简化为对式（2）求一阶导数并考虑到式（1）得挠曲线近似微分方程对式（4）求一阶导数并考虑到式（2）得对式（5）求一阶导数并考虑到式（1）得令，上式变为求解式（6）得挠曲线方程对式（7）求一阶导数得转角w′（x）满足方程1.2 变形边界条件在式（7）～式（8）中令x=0，得到A端的变形满足在式（7）～式（8）中令x=l，得到B端的变形满足1.3 内力边界条件对式（8）求一阶导数并考虑到式（4），得对式（13）求一阶导数并考虑到式（2），得在式（13）～式（14）中令x=0，得到A端的内力满足在式（13）～式（14）中令x=l，得到B端的内力满足1.4 特征方程根据表1所述不同支承情况的压杆，由变形和内力为零的边界条件，选用式（9）～式（12），式（15）～式（18）中的4个，构成关于初参数C1，C2，C3，C4的齐次线性方程组.因初参数C1，C2，C3，C4不全为0，则方程组的系数行列式必为0，由此得到临界载荷q cr的特征方程.1.4.1 固支--固支由式（9）～式（12）及表1所述该支承情况下为零的边界条件，得特征方程为1.4.2 固支--铰支由式（9）～式（11）和式（17）及表1所述该支承情况下为零的边界条件，得特征方程为1.4.3 固支--定向定向支承是指仅对压杆的转角作刚性约束，即该处压杆的转角为0，而不限制其挠度，即该处压杆的剪力为0.由式（9），式（10），式（12），式（18）及表1所述该支承情况下为0的边界条件，得特征方程为1.4.4 固支--自由由式（9），式（10），式（17），式（18）及表1所述该支承情况下为零的边界条件，得特征方程为1.4.5 铰支--铰支由式（9），式（11），式（15），式（17）及表1所述该支承情况下为零的边界条件，得特征方程为1.4.6 铰支--定向由式（9），式（12），式（15），式（18）及表1所述该支承情况下为零的边界条件，得特征方程为1.5 临界载荷近似解对1.4.1～1.4.6所述的特征?方程求解，若其有最小正数解（kl）cr，则由k3得临界载荷qcr=用软件如《Maple》等绘制1.4.1～1.4.6所述的6种压杆的临界载荷特征方程曲线，见图3，每条曲线在k l轴上的所有截距均为相应的特征方程的解，取第一截距即最小正数值（kl）cr计算临界载荷.软件分析表明，易失稳的压杆为图中虚线表示的情况.（1）1.4.3所述固支--定向支承压杆第一截距（kl）cr≈3.02，则该式表明，将大小为q cr l且沿压杆轴向均布施加的载荷，换算成大小为0.35833q cr l的集中载荷施将式（19）换算为欧拉公式形式，即加在压杆端部时，二者对压杆的稳定性影响是等效的.（2）1.4.4所述固支--自由支承和1.4.6所述“铰支--定向”支承压杆二者特征方程曲线重合，第一截距（k l）cr≈1.85，则将式（21）换算为欧拉公式形式，即该式表明，将大小为q cr l且沿压杆轴向均布施加的载荷，换算成大小为0.38969q cr l的集中载荷施加在压杆端部时，二者对压杆的稳定性影响是等效的. 设压杆处于临界状态时，轴向均布载荷为q，轴向集中载荷为F；设F cr为未考虑轴向均布载荷q时，按欧拉公式计算的临界载荷.不妨令F= K cr F cr，且称K cr 为F cr的修正因子.（1）对1.4.3所述固支--定向支承压杆采用载荷换算与叠加的方法，结合式（20）得q和F同时作用时，临界载荷经验公式为即，得其中（2）对1.4.4所述固支--自由支承和1.4.6所述铰支--定向支承压杆采用载荷换算与叠加的方法，结合式（22）得q和F同时作用时，临界载荷经验公式为即得其中由于文献［1-3］及Bessel函数法（源于文献［1］）仅研究了本文所述固支--自由支承的压杆，且研究范围是ql/F cr介于0与0.5之间，故这里仅就该可比情况进行比较.先分别按文献［1-3］、Bessel函数法及本文方法得到的公式或数据计算出F cr 的修正因子K cr（本文应按式（26）计算K cr），再计算出K cr的平均值cr，最后计算出各K cr与cr的偏差δcr（%），结果见表2.表2中带下划线及粗体数字表明：（1）当ql/F cr=0时，5种方法研究结果一致；（2）当ql/F cr在0.05左右时，文献［2-3］的研究结果与平均结果较吻合，文献［1］次之；（3）当ql/F cr介于0.10～0.20时，本文研究结果与平均结果较吻合，Bessel函数法次之；（4）当ql/F cr介于0.25～0.5时，本文研究结果与平均结果较吻合，文献［3］次之.由文献［1］知，工程中一般0.1≤ql/F cr≤0.25.由表2知，当0.1≤ql/F cr≤0.5时，本文研究结果与平均结果较吻合.对固支--定向、固支--自由和铰支--定向支承的压杆，q对其稳定性有明显影响.压杆处于临界状态时，q和F依次满足经验式（23）和式（25）.当F计算结果为负时，表示需要施加拉力，才可能避免压杆失稳.其他支承的压杆，q对其稳定性无明显影响.显然，对于完全理想支承的压杆，参照经验公式计算临界载荷是比较安全的，其中u为压杆的长度因数，其值取决于压杆的支承情况.【相关文献】1吴晓，杨立军，郑长成.轴向压力和均布载荷作用下压杆的屈曲研究.四川建筑科学研究，2008，34（5）：34-35，652张平占.BHA受压失稳的模型及临界钻压公式.煤田地质与勘探，1995，23（5）：59-623马肖，韩春民.具有较大初始位移的拼组柱在分布与集中荷载共同作用下的承载能力分析.四川建筑科学研究，2009，35（3）：27-29。

重心插值配点法求解分数阶Fredholm积分方程虎晓燕;韩惠丽【摘要】Barycentric interpolation collocation method was the combination and promotion of interpolation method and the collocation method.It was a new kind of numerical calculation method which had excellent numerical stability, high precision,and computational efficiency.This new method was used to solve fractional order Fredholm integral equation, and was deduced discrete calculate formula of the barycentric interpolation collocation method of fractional order Fredholm integral equation.Through theoretical analysis, it derived its existence and uniqueness of solutions and error analysis.Finally,some numerical examples were used to contrast equidistant nodes and the second Chebyshev nodes to demonstrate the effectiveness and precision of this method, and then the conditions that affect the precision was obtained.%重心插值配点法是插值法和配点法的结合和推广,它具有稳定性好、高精度和计算效率高等优点.主要运用高精度无网格重心插值配点法求解分数阶Fredholm积分方程.首先推导了基于分数阶Fredholm积分方程重心插值配点法的离散公式,然后通过理论分析得出其解的存在唯一性与误差分析,最后利用数值算例通过对等距节点与第二类Chebyshev节点的对比,验证了所用方法的高精度和可靠性,并得出影响精度的条件.【期刊名称】《郑州大学学报（理学版）》【年(卷),期】2017(049)001【总页数】7页(P17-23)【关键词】重心插值配点法;高精度;分数阶积分方程;无网格【作者】虎晓燕;韩惠丽【作者单位】宁夏大学数学统计学院宁夏银川 750021;宁夏大学数学统计学院宁夏银川 750021【正文语种】中文【中图分类】O175.6分数阶微积分有着广阔的应用领域和前景，如生物传热、粘弹性力学、图像信号和风力发电等. 由于求分数阶微积分方程的解析解存在诸多困难，所以重点研究它的数值解. 目前微积分方程的数值解法较多，如有限差分法[1]、有限元法[2]、小波方法[3]、拉格朗日乘子法[4] 、变分迭代法[5]、同伦摄动法[6]等，也有与统计力学相结合的随机行走方法[7]等. 上述方法较多研究整数阶微积分方程，而分数阶微积分方程更复杂，精确数值模拟分数阶微积分方程的困难更大，现有文献中研究分数阶微积分方程数值解的理论较少，而且很少考虑计算精度或收敛速度等问题. 高精度无网格重心插值配点法是一种新型求解分数阶积分方程的数值方法[8]，该方法相对其他数值方法而言应用前景广阔，克服了一般插值方法的振荡性和向前稳定性，具有计算精度高，不用划分网格，程序容易实现等诸多优点.重心插值配点法多用于数值模拟积分微分方程和偏微分方程. 文献[9]运用线性重心有理插值配置法求解了一阶线性双曲型初值边值问题.文献[10]运用两种重心型插值求解了带初值问题和边值问题的二阶线性微分方程，并将重心插值法与拟谱法以及微分求积法做了比较，得出重心插值法比其他数值方法精度更高.文献[11]使用直接线性重心有理插值(DRQ)求解了线性Fredholm积分函数，用间接线性重心有理插值法(IRQ)求解了线性Volterra积分函数，并给出了误差估计，从而说明了重心插值法的高效性.与以往插值配点法不同，本文将运用重心插值配点法求解分数阶Fredholm积分方程，从而推导出分数阶Fredholm积分方程重心插值配点法的离散计算公式. 通过证明解的存在唯一性、误差估计以及数值算例验证本文方法的高精度和可靠性. 设xj(j=0,1,…,n+1)是函数f(x)的n+1个不同的插值节点，其对应的函数值为fj.如果要使用多项式插值，则可以在多项式空间中求插值多项式p(x)，使得p(xj)=fj,j=0,1,…,n，而多项式空间的次数不超过n.定义1[12] 重心Lagrange插值公式由式(1)可知,分子和分母都包含重心权ωj,由式(2)可知,重心权只依赖于插值节点的分布，因此插值节点直接影响着重心权.本文主要采用等距节点和第二类Chebyshev节点来研究分数阶积分方程的数值解.定义2[13] 设函数f(t)为定义在J0=(0,+)上分段连续，且在J=[0,+)的任意有限子区间上可积. 对∀t>0、Re ν>0的任意复变量ν，函数f(t)的ν阶R-L分数阶定义,其中：ν>0；Γ(·)为Gamma函数.本文主要探讨如下形式的分数阶Fredholm积分方程，其中：f(x)为未知函数；g(x)∈C[a,b]为已知函数；K(x,t)为区间[a,b]×[a,b]上的、关于变量x,t的二元连续积分核函数.将积分方程的积分区间[a,b]离散为a=x0<x1<x2<…<xn=b，则令未知函数f(x)在离散节点x0,x1,x2,…,xn上的函数值为f0,f1,f2,…,fn，记fi=f(xi),i=0,1,2,…,n.用重心型插值来近似未知函数f(x)，式中，n.Lj(x)称为重心型插值基函数，ωj称为重心插值权，它只依赖于插值节点的分布. 将式(4)代入方程(3)得令方程(5)在节点x0,x1,x2,…,xn上均成立，可得令,则其中，且，从而式(6)可以写作将式(7)写成矩阵形式为其中:K=[Kij]:=[Kj(xi)]为未知函数f(x)在节点x0,x1,x2,…,xn的积分矩阵；I为n阶单位矩阵；F=[f(x0),f(x1),f(x2),…,f(xn)]T为f(x)在插值点处函数值；G=[g(x0),g(x1),g(x2),…,g(xn)]T为自由项g(x)所对应的插值点处的函数值. 这样式(8)就给出了积分方程(3)的重心插值配点法的离散计算公式.由Gauss公式[14]，对给定的xi有其中:tk,Ak分别为Gauss求积公式的积分点和积分权；q为积分点的数量，2为定积分的变换系数. Gauss求积公式是一种机械型求积方法，这类方法将积分求值问题转化为函数值的计算问题，积分权Ak紧密依赖于积分点tk的选取，与被积函数f(x)的形式没有关系，而积分点tk的数量q可以直接影响数值结果的好坏，当q越大计算误差越小，当q达到一定数量时，继续增加q对计算误差的影响将变小，但计算时间会变大，故可根据实际问题的需要适当选取q.4.1 解的存在唯一性引理1[15] (Fredholm选择定理)第二类Fredholm积分方程，要么对于所有函数f(x)都有解且解唯一，要么相应的齐次方程至少有一个非平凡解，即至少有一个不恒等于零的解.引理2[16] (压缩映射原理)设(X,d)是完备的距离空间， T：X→X是压缩映射.则T 在X中存在唯一的不动点x0.由引理1可知，本文所研究的方程必然有解. 在Banach空间中，设核K(x,t)在矩形区域D1=[a,b]×[a,b]上连续，自由项g(x)∈C[a,b]，下面利用压缩映射原理(引理2)证明其解的唯一性.定理1 方程(3)存在唯一解，当且仅当，其中t.证明定义映射T为则T:C[a,b]→C[a,b],由于x-t>0,从而(x-t)α-1>0,且对于任意的f1,f2∈C[a,b]，有‖Tf1-Tf2‖‖f1-f2‖‖f1-f2‖,式中令dt ，从而当<1时，T是C[a,b]上的压缩映射，容易得到C[a,b]是一个完备空间，从而可以利用引理2得知方程(3)存在唯一解，得证.4.2 误差估计本文主要采用两种范数误差来刻画数值解与精确解的逼近程度问题，其中L误差定义为‖ek‖，L2误差定义为.下面主要运用L误差研究本文方法的数值逼近问题.定理2 设(x)是方程(3)的近似解，若g(x)∈C[a,b]，且f(x)在[a,b]上存在n+1阶导数，则对∀x∈[a,b],有估计‖f(x)‖，式中，证明设(x)是方程(3)的近似解，则可得‖‖=‖‖≤ ,重心Lagrange插值公式 (1)可以写为p(x).若f(x)在[a,b]上存在n+1阶导数，则对∀x∈[a,b]，其重心插值有估计[17]其中:[a,b].现将式(10)代入式(9) ，可得方程(3)的解的误差估计‖‖.证毕.为了验证本文方法的高精度和可靠性，对2个具有精确解的数值算例进行数值模拟.算例中均分别选取等距节点和第二类Chebyshev节点，并对这两种节点下数值结果的L误差、L2误差和相对误差做比较. 下述算例的程序均在Matlab7.0中实现.算例1当α=1，取g(x)=11x/(12-x2)时,该分数阶方程退化为整数阶积分方程，精确解为f(x)=x(1-x). 表1为本文方法当α=1，n=10，q=8时两类节点和文献[18]方法的L误差、全局误差和相对误差对比，图1为α=1，n=10，q=8时两类节点误差与数值解、精确解比较. 当α=5/2时，适当选取g(x)使算例1精确解为f(x)=x(1-x).表2为α=5/2,n=10，q取不同值时两类节点的数值解、精确解、L误差、L2误差和相对误差. 其中全局误差定义为‖e‖2/n.由表1、表2和图1可知，用重心插值配点法求解分数阶Fredholm积分方程时,α增大，q增大时，其计算精度随之增大. 当q由10变为100时，精度增加了3个数量级，当q由100变为500时，精度增加了2个数量级，这说明q影响数值解的计算精度为,q越大计算效果越好，但相应计算时间增大. 两类节点所得数值解的计算精度都比较高，其中第二类Chebyshev节点比等距节点的计算精度更高. 当α=1时，算例1退化为整数阶时，q取10时相对误差达到10-16，与文献[18]的方法相比，本文方法的计算精度更高.算例2当α=5/2时，选取适当的g(x)使算例2的精确解为f(x)=x3(1-x)ex. 具体计算中，分别取11个和21个插值节点，即分为n=10和n=20两种情况，q取10. 由表3可以看出，Gauss积分点数量q相同时，n=20比n=10计算效果更好，等距节点比第二类Chebyshev节点更好.综合上述2个算例可知，应用重心插值配点法求解分数阶Fredholm积分方程时，可以得到很高的计算精度.计算精度随着阶数α的增大而增大，并且也受Gauss积分点数量q与插值节点个数n的影响. 当α=1，即积分方程由分数阶退化为整数阶，q与n取很小时就能得到很高的计算精度. 当α≠1，α越大，q越大，n越大时，分数阶积分方程的计算精度越高. 同时也可知，对于分数阶积分方程的精确解为基本初等函数时，第二类Chebyshev 节点比等距节点的精度更高，对于其精确解为更加复杂的函数时，等距节点比第二类Chebyshev 节点的精度更高一些.本文应用重心插值配点法求解了分数阶第二类Fredholm积分方程，推导出了分数阶Fredholm积分方程重心插值配点法的离散计算公式，证明了本文所讨论方程的解的存在唯一性，得出了误差估计，最后利用数值算例通过对比等距节点与第二类Chebyshev节点的数值结果，验证了本文方法的高精度和可靠性，并得出影响精度的条件，也给出两类节点的适用范围. 本文利用重心插值配点法对分数阶Fredholm积分方程进行了研究，得到了较好的结果.【相关文献】[1] 丁恒飞. 分数阶偏微分方程的有限差分方法[D]. 上海：上海大学，2014.[2] AGRAWAL O P. A general finite element formulation for fractional variational problems[J]. Journal of mathematical analysis and applications，2008，337(1)： 1-12. [3] HEYDARIA M H，HOOSHMANDASL M R，MOHAMMADI F，et al. Wavelets methodfor solving systems of nonlinear singular fractional Volterra integro-differential equations[J]. Communication in nonlinear science and numerical simulation，2014，19(1)：37-48.[4] 吴晓，黄志刚，杨立军. 用拉格朗日乘子法求解双模量静不定结构[J]. 力学与实践，2013，35(6)：79-81.[5] HETMANIOK E，NOWAK I，STOA D,et al. A study of the convergence of and error estimation for the homotopy perturbation method for the Volterra-Fredholm integral equations[J]. Applied mathematics letters，2013，26(26)：165-169.[6] ZHAN J，CHEN L，FU Y, et al． A new random walk simulation model for study of diffusion behavior of single particle within two-dimensional space[J]. Journal of electrochemistry，2012，18(5)：427-436.[7] WU G C. A fractional variational iteration method for solving fractional nonlinear differential equations[J]. Computers and mathematics with applications，2011，61(61)：2186-2190.[8] 王兆清，李淑萍，唐炳涛.一维重心型插值：公式、算法和应用[J].山东建筑大学学报，2007，22(5)：448-453.[9] BALTENSPERGER R, BERRUT J P. The linear rational collocation method[J]. Journal of computational and applied mathematics，2001，134(1/2)：243-258.[10]李淑萍. 基于重心型插值的数值计算方法[J]. 山东科学，2010，23(4)：13-16.[11]KLEIN G，BERRUT J P. Linear barycentric rational quadrature[J]. Bit numerical mathematics，2012，52(2)：407-424.[12]李树忱，王兆清. 高精度无网格重心插值配点法:算法、程序及工程应用[M]. 北京：科学出版社，2012.[13]OLDHAM K B，SPANIER J. The fractional calculus theory and applications of differentiation and integration to arbitrary order[M]. New York：Dover Publications Inc，2002.[14]李岳生，黄友谦. 数值逼近[M]. 北京：人民教育出版社，1978.[15]吕涛，黄晋. 积分方程的高精度算法[M]. 北京：科学出版社，2013.[16]匡继昌. 实分析与泛函分析[M]. 北京：高等教育出版社，2002.[17]ZOU L,LI C. Barycentric Lagrange blending rational interpolation based on Padè approximation[C]//Proceeding of the 2011 International Conference on Computational and Information Sciences.Washington,2011：1124-1127.[18]张倩，韩惠丽，张盼盼. 基于有理 Haar小波求解分数阶第2类Fredholm积分方程 [J]. 江西师范大学学报(自然科学报)，2014，38(1)：47-50.。

轴向均布荷载压杆稳定问题的重心插值配点法
于卫涛;宋洁;赵维霞
【期刊名称】《山东建筑大学学报》
【年(卷),期】2011(026)004
【摘要】用重心插值配点法对轴向均布荷载下压杆稳定问题进行了研究.采用重心Lagrange插值多项式建立未知函数的微分矩阵,采用配点法将压杆稳定问题的控制方程表示为代数方程组.通过求解代数方程组系数矩阵的特征值和特征向量,得到了精度很高的后屈曲挠度数值和临界载荷数值.实例证明,这种方法原理简单,易于程序实现.
【总页数】4页(P353-355,360)
【作者】于卫涛;宋洁;赵维霞
【作者单位】滨州市建设工程质量监督站,山东滨州256600;滨州市工程建设标准定额站,山东滨州256600;山东城市建设职业学院,山东济南250014
【正文语种】中文
【中图分类】TU311.41
【相关文献】
1.轴向均布载荷下压杆稳定问题的DQ解 [J], 刘洋;杨永波;梁枢平
2.重心有理插值配点法分析压杆稳定问题 [J], 王奇;王兆清
3.重心Lagrange插值配点法求解二维非线性Volterra积分方程 [J], 张薇
4.重心插值配点法求解二维Sobolev方程 [J], 宋灵宇;武莉莉;卢梦双
5.重心插值配点法求解Cahn-Hilliard方程 [J], 邓杨芳;黄蓉;翁智峰
因版权原因，仅展示原文概要，查看原文内容请购买。