论文部分内容阅读
摘 要：工程领域中随机振动是一类非常普遍的现象，构件动力响应常需要计算一类周期函数的广义积分，目前现有的方法表达式复杂，计算精度和效率低下.利用周期函数的特点将广义积分转化为有限区间积分，获得简明的封闭解，然后根据高斯-切比雪夫积分具有计算精度和效率高的特点，推导出该类广义积分的新近似解.通过算例对比分析，验证了本文所提方法的正确性和高效性，对解决工程领域的振动响应分析具有重要的参考价值.
关键词：周期函数 ;广义积分; 高斯-切比雪夫积分; 简明近似解
中图分类号：TU311.3;O21 DOI：10.16375/j.cnki.cn45-1395/t.2021.03.001
0 引言
随机振动存在于各类工程中，是工程设计中构件安全设计必须考虑的重要因素[1-3]，而频域法是求解各类工程结构随机动力响应的重要方法[4-5].李创第等[6]利用虚拟激励法研究了非粘滞阻尼结构基于随机激励下结构响应的谱矩，需要对结构响应功率谱在[0，+∞）的积分区间进行数值积分.李暾等[7]利用二次正交化法研究了建筑结构基于近似Davenport风速谱的结构随机风振响应的简明封闭解，有效地提高了随机响应分析的精度和效率.而分析非平稳随机激励[2-3]或考虑行波效应的平稳激励[8]下的结构随机响应时需要计算一类含有周期函数的2个积分：
[0∞sinωtp2+ω2dω] （1）
[0∞ωcosωtp2+ω2dω] （2）
其中：[t≥0]且為实数，p为常数且实数或者复数.
目前对于上述2个积分的解法有基于留数定律[9-10]的解法，该解法异常复杂，将上述积分表示成指数积分的封闭解，且只能用于p为实常数的情况，而指数积分仍然为数值解.此外，针对上述2个积分，Newton-Cotes积分是一种直接方法[11]，其积分精度受积分步长和积分上限的影响较大，且目前未见有积分上限取值的研究文献.高斯系列积分[11]利用带权的正交多项式来计算积分，具有计算精度和计算效率高的优点，有着广泛的应用.目前利用高斯积分研究式（1）、式（2）2个函数的积分的文献鲜有.
针对式（1）、式（2）2个函数的积分，首先利用其周期性和收敛性的特点，将广义积分转变为有限个[0，[2π]]的积分，然后利用三角函数的特点，将 [0，[2π]]的积分变化为高斯-切比雪夫积分，从而提出了一种计算2个积分的新简明解.
1 简明近似解
1.1 式（1）的近似解
利用sin函数的周期性，作如下变换：
[0∞sinωtp2+ω2dω=x=ωtt0∞sinxt2p2+x2dx=]
[tk=0∞02πsinxt2p2+x+2πk2dx] （3）
由式（1）可知，随着k的增加，积分项的值收敛于0，因此，实际应用时k值取有限值.同时，从式（3）可知，积分上限由正无穷大换成[2π]，利用sin函数的特点，式（3）改写为：
[0∞sinωtp2+ω2dω=tk=0∞Ak+Bk-tC] （4）
式中：
[Ak=-π2π2 sinxt2p2+x+2πk2dx=y=sinx-11yt2p2+arcsiny+2πk21-y2dyBk=-π2π2 sinxt2p2+x+2πk+π2dx=y=sinx-11 yt2p2+arcsiny+2πk+π21-y2dy]（5）
[C=-π20 sinxt2p2+x2dx] （6）
式（5）中的[Ak]及[Bk]可统一表示为：
[Ak（Bk）=-11 （-1）kyt2p2+arcsiny+πk21-y2dy]（7）
高斯系列数值积分[11]具有计算速度快、精度高的特点，而式（7）满足高斯-切比雪夫积分的积分区间[-1， 1]和权值函数[1-y2-0.5]，则式（7）的高斯积分为[8]：
[Ak（Bk）=-11 （-1）kyt2p2+arcsiny+πk21-y2dy=]
[i=1n（-1）kyit2p2+arcsinyi+πk2ai] （8）
式中：[ai]为高斯求积系数，[i]为高斯积分点数.
对式（6）采用复化Simpson公式计算：
[C=-π20sinxt2p2+x2dx=]
[h3-1t2p2+π2/4+4i=1msinx2it2p2+x22i+2i=1msinx2i-1t2p2+x22i-1] （9）
式中：[h=π4m]，[xi=-π2+ih]，[m]为节点数.
把式（8）、式（9）代入式（4），得出： [0∞sinωtp2+ω2dω=tk=0∞i=1n（-1）kyit2p2+arcsinyi+πk2-th3-1t2p2+π2/4+4i=1msinx2it2p2+x22i+2i=1msinx2i-1t2p2+x22i-1] （10）
1.2 式（2）的近似解
利用cos函数的周期性，作如下变换：
[0∞ωcosωtω2+p2dω=x=ωt0∞xcosxx2+tp2dx=]
[k=0∞02π（2kπ+x）cos（x+2kπ）（2kπ+x）2+tp2dx] （11）
由式（11）可知，随着k的增加，积分项的值收敛于0，因此，实际应用时k值取有限值.同时，从式（11）可知，积分上限可由正无穷大换成2p.利用cos函数的特点，式（11）改写为：
[0∞ωcosωtω2+p2dω=k=0∞0π（2kπ+x）cosx（2kπ+x）2+tp2dx+]
[k=0∞π2π（2kπ+x）cosx（2kπ+x）2+tp2dx] （12）
利用cos函数的特点，式（12）进一步改写为：
[0∞ωcosωtω2+p2dω=k=0∞0π（2kπ+x）cosx（2kπ+x）2+tp2dx-]
[k=0∞0π（2kπ+x+π）cosx（2kπ+x+π）2+tp2dx] （13）
最后式（13）简化为：
[0∞ωcosωtω2+p2dω=]
[k=0∞（-1）k+1π0（kπ+x）cosx（kπ+x）2+tp2dx] （14）
对式（14）积分变换：
[0∞ωcosωtω2+p2dω=]
[k=0∞（-1）k+1-11-（kπ+arccosy）y（kπ+arccosy）2+tp21-y2dy] （15）
式（15）满足高斯-切比雪夫积分的积分区间 [-1， 1]，且权值函数为[1-y2-0.5]，则式（15）的高斯积分为[11]：
[0∞ωcosωtω2+p2dω=]
[k=0∞（-1）ki=1n（kπ+arccosyi）yiai（kπ+arccosyi）2+tp2] （16）
式中：[ai]为高斯求积系数，[i]为高斯积分点数.
2 数值实验
为验证本文方法的计算精度和效率，以 [p=-3+4i （i=-1）]为例，利用基于常规的数值积分方法和本文方法进行计算和说明.
[0∞sin4ω-3+4i2+ω2dω] （17）
[0∞ωcos4ω-3+4i2+ω2dω] （18）
對式（17）和式（18）采用梯形积分时，其表达式可表示为：
[0∞sinωtp2+ω2dω=]
[i=0n（12（sinωitp2+ω2i+sinωi+1tp2+ω2i+1）Δω）] （19）
[0∞ωcosωtp2+ω2dω=]
[i=0n（12（ωicosωitp2+ω2i+ωi+1cosωi+1tp2+ω2i+1）Δω）] （20）
式中：[ωi=iΔω]，[Δω]为积分点间距，i为整数;n为积分上限值参数.式（19）及式（20）的计算精度和效率受积分上限和积分间距2个因素影响，而对于本文方法则受积分上限和高斯积分点个数的影响较大.为此，就上述2个因素进行分析，具体见表1—表4.
从表1来看，式（17）利用本文方法计算时，当上限[z≥]50，节点数[≥]8，所获得积分结果已趋于稳定.从表2可知，传统数值积分法的积分上限[k≥]500 rad/s，积分间距[≤]0.001 0 rad/s，积分结果也趋于稳定.比较稳定后的数值可知，本文方法与传统数值积分法非常接近，但本文方法在 [z=]50，节点数=8时，耗时为0.112 s;而传统积分方法耗时在[k=]500 rad/s，积分间距=0.001 0 rad/s，耗时为0.315 s，说明了本文方法的效率比传统积分方法更高.此外，传统积分方法的积分上限和积分间距需要在较大范围内去试算才能确定，而本文方法节点数和k值均较小，说明本文方法计算稳定性好.
从表3来看，式（18）利用本文方法计算时，当上限[z≥]200，节点数[≥]10所获得积分结果已趋于稳定.从表4可知，传统数值积分法的积分上限 [k≥]50 000 rad/s，积分间距[≤]0.001 0 rad/s，积分结果也趋于稳定.比较稳定后的数值可知，本文方法与传统数值积分法非常接近，但本文方法在[z=]200，节点数=10时，耗时为0.407 s;而传统积分方法在[k=]50 000 rad/s，积分间距=0.010 0 rad/s时，耗时为3.887 s，说明了本文方法的效率比传统积分方法更高.此外，传统积分方法的积分上限和积分间距需要在较大范围内去试算才能确定，而本文方法节点数和k值均较小，说明本文方法计算稳定性好. 3 结论
本文针对随机振动中两类广义积分[0∞sinωtp2+ω2dω]和[0∞ωcosωtp2+ω2dω]无简明解的问题，利用周期函数和高斯-切比雪夫积分提出了一种计算精度和效率高的简明近似解.此外，传统积分方法的积分上限和积分间距需要在较大范围内去试算才能确定，而本文方法节点数和k值均在较小的范围内试算就能确定，说明本文方法计算上述两类广义积分具有较好的稳定性.因此，本文方法对解决工程领域的振动响应分析具有重要的参考价值.
参考文献
[1] CRANDALL S H，MARK W D. Random vibration in mechanical systems[M].New York：Academic Press，1963.
[2] BARBATO M，CONTE J P. Time-variant reliability analysis of linear elastic systems subjected to fully nonstationary stochastic excitations[J].Journal of Engineering Mechanics，2015，141（S6）：1-10.
[3] PENG B F，CONTE J P. Closed-form solutions for the response of linear systems to fully nonstationary earthquake excitation[J]. Journal of Engineering Mechanics，1998，124（6）：684-694.
[4] 郑兆昌.随机振动矩阵直接谱分析法[C]//第二十三届全国振动与噪声控制学术会议， 沈阳，2010.
[5] 林家浩，张亚辉，赵岩.虚拟激励法在国内外工程界的应用回顾与展望[J].应用数学和力学，2017，38（1）：1-32.
[6] 李創第，贺王涛，葛新广.卷积型非粘滞阻尼结构随机地震动系列响应求解的虚拟激励法[J].广西科技大学学报，2021，32（1）：78-84.
[7] 李暾，张梦丹，姜琰，等.基于近似Davenport风速谱的建筑结构动力响应的新封闭解法[J].广西科技大学学报，2020，31（4）：1-10，18.
[8] 张文首，林家浩.大跨度结构考虑行波效应时平稳随机地震响应的闭合解[J].固体力学学报，2004（4）：446-450.
[9] 王春.留数定理在Mellin逆变换中的应用[J].大学数学，2020，36（2）：106-110.
[10] 数学手册编写组.数学手册[M].北京：人民教育出版社，1979.
[11] 朱建新，李有法. 数值计算方法[M].北京：高等教育出版社，2012.
Numerical method for generalized integral of a class of
periodic functions
GE Xinguang， LI Yuxiang， YANG Xuefeng
（School of Civil Engineering and Architecture， Guangxi University of Science and Technology，
Liuzhou 545006， China）
Abstract： Random vibration is a very common phenomenon in the engineering field. The generalized integral of a class of periodic functions is often needed to calculate the dynamic response of components. The existing methods have complex expressions and low calculation accuracy and efficiency. In this paper， the generalized integral is transformed into a finite interval integral by using the characteristics of periodic function， and a concise closed solution is obtained. Then， a new approximate solution of this kind of integral is derived by combining the high accuracy and efficiency of Gauss-Chebyshev integral. The correctness and efficiency of the proposed method are verified through the comparative analysis of examples， which has important value for solving the vibration response analysis in engineering field.
Key words： periodic function; generalized integral; Gauss-Chebyshev integral; concise approximate solution
（责任编辑：罗小芬、黎 娅）
关键词：周期函数 ;广义积分; 高斯-切比雪夫积分; 简明近似解
中图分类号：TU311.3;O21 DOI：10.16375/j.cnki.cn45-1395/t.2021.03.001
0 引言
随机振动存在于各类工程中，是工程设计中构件安全设计必须考虑的重要因素[1-3]，而频域法是求解各类工程结构随机动力响应的重要方法[4-5].李创第等[6]利用虚拟激励法研究了非粘滞阻尼结构基于随机激励下结构响应的谱矩，需要对结构响应功率谱在[0，+∞）的积分区间进行数值积分.李暾等[7]利用二次正交化法研究了建筑结构基于近似Davenport风速谱的结构随机风振响应的简明封闭解，有效地提高了随机响应分析的精度和效率.而分析非平稳随机激励[2-3]或考虑行波效应的平稳激励[8]下的结构随机响应时需要计算一类含有周期函数的2个积分：
[0∞sinωtp2+ω2dω] （1）
[0∞ωcosωtp2+ω2dω] （2）
其中：[t≥0]且為实数，p为常数且实数或者复数.
目前对于上述2个积分的解法有基于留数定律[9-10]的解法，该解法异常复杂，将上述积分表示成指数积分的封闭解，且只能用于p为实常数的情况，而指数积分仍然为数值解.此外，针对上述2个积分，Newton-Cotes积分是一种直接方法[11]，其积分精度受积分步长和积分上限的影响较大，且目前未见有积分上限取值的研究文献.高斯系列积分[11]利用带权的正交多项式来计算积分，具有计算精度和计算效率高的优点，有着广泛的应用.目前利用高斯积分研究式（1）、式（2）2个函数的积分的文献鲜有.
针对式（1）、式（2）2个函数的积分，首先利用其周期性和收敛性的特点，将广义积分转变为有限个[0，[2π]]的积分，然后利用三角函数的特点，将 [0，[2π]]的积分变化为高斯-切比雪夫积分，从而提出了一种计算2个积分的新简明解.
1 简明近似解
1.1 式（1）的近似解
利用sin函数的周期性，作如下变换：
[0∞sinωtp2+ω2dω=x=ωtt0∞sinxt2p2+x2dx=]
[tk=0∞02πsinxt2p2+x+2πk2dx] （3）
由式（1）可知，随着k的增加，积分项的值收敛于0，因此，实际应用时k值取有限值.同时，从式（3）可知，积分上限由正无穷大换成[2π]，利用sin函数的特点，式（3）改写为：
[0∞sinωtp2+ω2dω=tk=0∞Ak+Bk-tC] （4）
式中：
[Ak=-π2π2 sinxt2p2+x+2πk2dx=y=sinx-11yt2p2+arcsiny+2πk21-y2dyBk=-π2π2 sinxt2p2+x+2πk+π2dx=y=sinx-11 yt2p2+arcsiny+2πk+π21-y2dy]（5）
[C=-π20 sinxt2p2+x2dx] （6）
式（5）中的[Ak]及[Bk]可统一表示为：
[Ak（Bk）=-11 （-1）kyt2p2+arcsiny+πk21-y2dy]（7）
高斯系列数值积分[11]具有计算速度快、精度高的特点，而式（7）满足高斯-切比雪夫积分的积分区间[-1， 1]和权值函数[1-y2-0.5]，则式（7）的高斯积分为[8]：
[Ak（Bk）=-11 （-1）kyt2p2+arcsiny+πk21-y2dy=]
[i=1n（-1）kyit2p2+arcsinyi+πk2ai] （8）
式中：[ai]为高斯求积系数，[i]为高斯积分点数.
对式（6）采用复化Simpson公式计算：
[C=-π20sinxt2p2+x2dx=]
[h3-1t2p2+π2/4+4i=1msinx2it2p2+x22i+2i=1msinx2i-1t2p2+x22i-1] （9）
式中：[h=π4m]，[xi=-π2+ih]，[m]为节点数.
把式（8）、式（9）代入式（4），得出： [0∞sinωtp2+ω2dω=tk=0∞i=1n（-1）kyit2p2+arcsinyi+πk2-th3-1t2p2+π2/4+4i=1msinx2it2p2+x22i+2i=1msinx2i-1t2p2+x22i-1] （10）
1.2 式（2）的近似解
利用cos函数的周期性，作如下变换：
[0∞ωcosωtω2+p2dω=x=ωt0∞xcosxx2+tp2dx=]
[k=0∞02π（2kπ+x）cos（x+2kπ）（2kπ+x）2+tp2dx] （11）
由式（11）可知，随着k的增加，积分项的值收敛于0，因此，实际应用时k值取有限值.同时，从式（11）可知，积分上限可由正无穷大换成2p.利用cos函数的特点，式（11）改写为：
[0∞ωcosωtω2+p2dω=k=0∞0π（2kπ+x）cosx（2kπ+x）2+tp2dx+]
[k=0∞π2π（2kπ+x）cosx（2kπ+x）2+tp2dx] （12）
利用cos函数的特点，式（12）进一步改写为：
[0∞ωcosωtω2+p2dω=k=0∞0π（2kπ+x）cosx（2kπ+x）2+tp2dx-]
[k=0∞0π（2kπ+x+π）cosx（2kπ+x+π）2+tp2dx] （13）
最后式（13）简化为：
[0∞ωcosωtω2+p2dω=]
[k=0∞（-1）k+1π0（kπ+x）cosx（kπ+x）2+tp2dx] （14）
对式（14）积分变换：
[0∞ωcosωtω2+p2dω=]
[k=0∞（-1）k+1-11-（kπ+arccosy）y（kπ+arccosy）2+tp21-y2dy] （15）
式（15）满足高斯-切比雪夫积分的积分区间 [-1， 1]，且权值函数为[1-y2-0.5]，则式（15）的高斯积分为[11]：
[0∞ωcosωtω2+p2dω=]
[k=0∞（-1）ki=1n（kπ+arccosyi）yiai（kπ+arccosyi）2+tp2] （16）
式中：[ai]为高斯求积系数，[i]为高斯积分点数.
2 数值实验
为验证本文方法的计算精度和效率，以 [p=-3+4i （i=-1）]为例，利用基于常规的数值积分方法和本文方法进行计算和说明.
[0∞sin4ω-3+4i2+ω2dω] （17）
[0∞ωcos4ω-3+4i2+ω2dω] （18）
對式（17）和式（18）采用梯形积分时，其表达式可表示为：
[0∞sinωtp2+ω2dω=]
[i=0n（12（sinωitp2+ω2i+sinωi+1tp2+ω2i+1）Δω）] （19）
[0∞ωcosωtp2+ω2dω=]
[i=0n（12（ωicosωitp2+ω2i+ωi+1cosωi+1tp2+ω2i+1）Δω）] （20）
式中：[ωi=iΔω]，[Δω]为积分点间距，i为整数;n为积分上限值参数.式（19）及式（20）的计算精度和效率受积分上限和积分间距2个因素影响，而对于本文方法则受积分上限和高斯积分点个数的影响较大.为此，就上述2个因素进行分析，具体见表1—表4.
从表1来看，式（17）利用本文方法计算时，当上限[z≥]50，节点数[≥]8，所获得积分结果已趋于稳定.从表2可知，传统数值积分法的积分上限[k≥]500 rad/s，积分间距[≤]0.001 0 rad/s，积分结果也趋于稳定.比较稳定后的数值可知，本文方法与传统数值积分法非常接近，但本文方法在 [z=]50，节点数=8时，耗时为0.112 s;而传统积分方法耗时在[k=]500 rad/s，积分间距=0.001 0 rad/s，耗时为0.315 s，说明了本文方法的效率比传统积分方法更高.此外，传统积分方法的积分上限和积分间距需要在较大范围内去试算才能确定，而本文方法节点数和k值均较小，说明本文方法计算稳定性好.
从表3来看，式（18）利用本文方法计算时，当上限[z≥]200，节点数[≥]10所获得积分结果已趋于稳定.从表4可知，传统数值积分法的积分上限 [k≥]50 000 rad/s，积分间距[≤]0.001 0 rad/s，积分结果也趋于稳定.比较稳定后的数值可知，本文方法与传统数值积分法非常接近，但本文方法在[z=]200，节点数=10时，耗时为0.407 s;而传统积分方法在[k=]50 000 rad/s，积分间距=0.010 0 rad/s时，耗时为3.887 s，说明了本文方法的效率比传统积分方法更高.此外，传统积分方法的积分上限和积分间距需要在较大范围内去试算才能确定，而本文方法节点数和k值均较小，说明本文方法计算稳定性好. 3 结论
本文针对随机振动中两类广义积分[0∞sinωtp2+ω2dω]和[0∞ωcosωtp2+ω2dω]无简明解的问题，利用周期函数和高斯-切比雪夫积分提出了一种计算精度和效率高的简明近似解.此外，传统积分方法的积分上限和积分间距需要在较大范围内去试算才能确定，而本文方法节点数和k值均在较小的范围内试算就能确定，说明本文方法计算上述两类广义积分具有较好的稳定性.因此，本文方法对解决工程领域的振动响应分析具有重要的参考价值.
参考文献
[1] CRANDALL S H，MARK W D. Random vibration in mechanical systems[M].New York：Academic Press，1963.
[2] BARBATO M，CONTE J P. Time-variant reliability analysis of linear elastic systems subjected to fully nonstationary stochastic excitations[J].Journal of Engineering Mechanics，2015，141（S6）：1-10.
[3] PENG B F，CONTE J P. Closed-form solutions for the response of linear systems to fully nonstationary earthquake excitation[J]. Journal of Engineering Mechanics，1998，124（6）：684-694.
[4] 郑兆昌.随机振动矩阵直接谱分析法[C]//第二十三届全国振动与噪声控制学术会议， 沈阳，2010.
[5] 林家浩，张亚辉，赵岩.虚拟激励法在国内外工程界的应用回顾与展望[J].应用数学和力学，2017，38（1）：1-32.
[6] 李創第，贺王涛，葛新广.卷积型非粘滞阻尼结构随机地震动系列响应求解的虚拟激励法[J].广西科技大学学报，2021，32（1）：78-84.
[7] 李暾，张梦丹，姜琰，等.基于近似Davenport风速谱的建筑结构动力响应的新封闭解法[J].广西科技大学学报，2020，31（4）：1-10，18.
[8] 张文首，林家浩.大跨度结构考虑行波效应时平稳随机地震响应的闭合解[J].固体力学学报，2004（4）：446-450.
[9] 王春.留数定理在Mellin逆变换中的应用[J].大学数学，2020，36（2）：106-110.
[10] 数学手册编写组.数学手册[M].北京：人民教育出版社，1979.
[11] 朱建新，李有法. 数值计算方法[M].北京：高等教育出版社，2012.
Numerical method for generalized integral of a class of
periodic functions
GE Xinguang， LI Yuxiang， YANG Xuefeng
（School of Civil Engineering and Architecture， Guangxi University of Science and Technology，
Liuzhou 545006， China）
Abstract： Random vibration is a very common phenomenon in the engineering field. The generalized integral of a class of periodic functions is often needed to calculate the dynamic response of components. The existing methods have complex expressions and low calculation accuracy and efficiency. In this paper， the generalized integral is transformed into a finite interval integral by using the characteristics of periodic function， and a concise closed solution is obtained. Then， a new approximate solution of this kind of integral is derived by combining the high accuracy and efficiency of Gauss-Chebyshev integral. The correctness and efficiency of the proposed method are verified through the comparative analysis of examples， which has important value for solving the vibration response analysis in engineering field.
Key words： periodic function; generalized integral; Gauss-Chebyshev integral; concise approximate solution
（责任编辑：罗小芬、黎 娅）