一类含符号函数的无穷积分的高效数值解法

葛新广 卢嘉康 张梨荣 罗臻

摘 要：基于频域法研究几类随机激励下工程结构的动力响应，需要求解4种含符号函数在无穷区间的积分解，而该4种积分式目前存在计算效率和精度的问题。首先，根据积分运算法则将4种积分式的计算转换为一种积分式的计算；其次，基于留数定理和高斯-雅克比积分提出了4种积分式的高效简明数值解法；最后，通过算例研究了传统方法受积分上限和积分区间的影响、本文方法受高斯-雅克比积分节点数的影响。结果表明：传统方法的计算结果受积分区间和积分上限的影响较大，而本文方法取25个高斯-雅克比积分节点数即可获得较高的精度，且具有较高的计算效率。

关键词：无穷积分；符号函数；高斯-雅克比积分；高效数值解

中图分类号：TU311.3 DOI：10.16375/j.cnki.cn45-1395/t.2024.03.007

0 引言

随机振动现象广泛存在于工程领域，如航空器在飞行过程中因风阻导致的振动[1]、土木工程领域的结构地震动[2]和风振动[3]、车辆在路面运行过程中因路面不平顺导致的振动[4-5]等。频域法是分析各类随机激励下结构构件动力响应的重要方法之一[6-7]。林家浩等[8]基于频域法提出了虚拟激励法，该方法广泛应用于各种结构的随机激励下的响应分析[9]，是研究随机地震动激励下土木工程结构的主要方法之一。实际地震动过程具有典型的非平稳特性，土木工程结构基于非平稳地震动激励下的动力响应分析[10]是评估结构地震过程安全的重要内容。大跨度土木工程结构在地震动激励下具有典型的行波效应[11]，考虑行波效应比不考虑行波效应的此类结构地震响应更加显著。火车基于轨道不平顺激励[12]和汽车基于路面不平顺激励下动力响应[13]分析需要考虑相干函数，考虑相干函数的影响会比不考虑行波效应的交通工具更加显著。上述工程领域的分析需要计算一类含有符号函数[sgn（ω）egiωtkiω-q-n]（sgn为符号函数，[ω]为积分变量，g=±1，i=[-1]，t为振动时长，k=±1，q为实部为负值的复数或者负实数，n为正整数）的无穷积分[14-16]，而目前针对这种积分所采用的方法较为复杂[15]，或者采用数值解的计算精度不易控制。张志超等[17]基于虚拟激励法-精细积分法研究了车桥耦合系统的非平稳激励下的车辆系统响应，所得解为数值解，其分析精度受时间步长和频率上限的影响较大，需要试算才能确定精度。Barbato等[14]研究了非平稳随机地震动激励下结构时变可靠度的非几何谱矩的计算，所得封闭解含有积分表达式，而该积分表达式是针对实部较小的复特征值，容易出现偏差较大的结果。葛新广[16]在Barbato[14]所提方法的基础上提出了上述4个表达式的改进计算方法，并将其应用于分数导数阻尼器耗能结构基于Conte-Peng完全非平稳激励下的响应分析，但所获得的表达式依然复杂，且存在对q值的实部较小时计算失真的情况。因此，有必要针对此类含有符号函数的积分提出更高效的计算方法。

高斯-雅克比积分法具有分析精度高和效率高的特点，广泛应用于各种复杂表达式的积分计算[16，18]。本文根据[sgn（ω）egiωtkiω-q-n]式中g、k的取值，首先研究了4种式子之间的关系式，并将4种积分式的计算转换为一种积分式的计算；其次，基于留数定理和高斯-雅克比积分提出了一种积分式的表达式的解，同时将其推广至其他3种表达式的计算；最后通过算例验证了所提方法的正确性和高效性。

1 4个无穷积分表达式的关系

根据g、k的取值，可将[sgn（ω）egiωtkiω-q-n]展开为4个无穷积分表达式

[A1q=-∞∞sgn（ω）eiωtiω-qndω ，] （1）

[A2q=-∞∞sgn（ω）e-iωtiω-qndω]， （2）

[A3q=-∞∞sgn（ω）eiωt-iω-qndω]， （3）

[A4q=-∞∞sgn（ω）e-iωt-iω-qndω]. （4）

对式（2）积分进行变换得

[A2q=-∞∞sgn（-ω）eiωt-iω-qnd-ω]. （5）

调换积分上下限及利用符号函数性质，式（5）改写为

[A2q=--∞∞sgn（ω）eiωt-iω-qndω]. （6）

为获得A2（q）与A1（q）的关系，对式（6）进一步改写为

[A2q=（-1）n-1-∞∞sgn（ω）eiωtiω+qndω]. （7）

对比式（7）与式（1）的关系，则存在

[A2q=（-1）n-1A1-q]. （8）

式（8）给出了A2（q）与A1（q）的关系式。

对式（3）进行积分变换可得

[A3q=-∞∞sgn（-ω）e-iωtiω-qnd-ω]. （9）

调整积分上下限并利用符号函数性质，对式（9）改写为

[A3q=--∞∞sgn（ω）e-iωtiω-qndω]. （10）

比较式（2）与式（10），则

[A3q=-A2q]. （11）

利用式（8），则建立A1（q）与A3（q）的关系式

[A3q=（-1）nA1-q]. （12）

对式（4）进行积分变换可得

[A4q=-∞∞sgn（-ω）eiωtiω-qnd-ω]. （13）

调整积分上下限并利用符号函数性质，对式（13）改写为

[A4q=--∞∞sgn（ω）eiωtiω-qndω]. （14）

比较式（14）与式（1），则

[A4q=-A1q]. （15）

式（15）给出了A1（q）与A4（q）的关系式。

式（8）、式（12）、式（15）、A2（q）—A4（q）均与A1（q）建立了关系，为此若获得式（1）—式（4）的积分解，只需要获得式（1）的积分解。

2 4个无穷积分的近似解

2.1 A1（q）的近似解

对式（1）进行改写得

[A1q=-in-∞∞sgn（ω）eiωtω+iqndω=-inA11+A12，]  （16）

式中：

[A11=0∞eiωtω+iqndω]， （17）

[A12=-∞0eiωtω+iqndω]. （18）

利用留数定理并结合式（17）及式（18）的特点，A11和A12的计算如图1所示。

[C5][C1][C3][C4][C2][实轴][虚轴] [O]

图1 留数定理应用示意图

[A11=0∞eiωtω+iqndω=C1eiωtω+iqndω，] （19）

[A12=-∞0eiωtω+iqndω=C5eiωtω+iqndω]. （20）

由留数定理，可构造积分函数为

[（21）]

[（22）]

[[式中：Rey=max0，ReyRey；                 Imy=max0，ImyImy；] ]

Re（y）、Im（y）分别表示求y的实部和虚部。

由式（21）及式（22），则A11与A12可表示为

[（23）]

[（24）]

对式（23）及式（24），C2—C4积分为

[C2eiztz+iqndz=limr→∞0π2eireiθtrieiθreiθ+iqndθ=0]，   （25）

[（26）]

[C4eiztz+iqndz=limr→∞π2πeireiθtrieiθreiθ+iqndθ=0]，    （27）

式中：r为留数定理中的半径。

把式（25）—式（27）代入式（23）和式（24），则A11和A12表示为

[A11=Re-iqIm-iq2πintn-1eqtn-1！+-in-1D]，   （28）

[A12=-2πintn-1eqtn-1！ReiqIm-iq+-in-1D]，   （29）

式中：

[D=0∞e-yty+qndy]. （30）

对式（30）利用高斯-雅克比积分进行高精度近似计算，首先对积分区间进行变换，令[x=1-y1+y]，则[y=1-x1+x]，代入式（30），则

[D=-11e-1-x1+xt1-x1+x+qn-21+x2dx]. （31）

对式（31）进一步简化可得

[D=2-11（1+x）n-2e-1-x1+xt1+q+q-1xndx]. （32）

利用高斯-雅克比积分对式（30）进行改写，即

[D≈2j=1N（1+xj）n-2Aje-1-xj1+xjt1+q+q-1xjn ，] （33）

式中：[xj、Aj]为权函数[ρ=1-x01+x0]的高斯-雅克比积分高斯点及系数；N为积分节点数。

把式（28）、式（29）及式（33）代入式（16），则

[A1q≈Im-iqRe-iq-Reiq2πtn-1eqtn-1！+4-1nij=1N（1+xj）n-2Aje-1-xj1+xjt1+q+q-1xjn.] （34）

由于工程上的q为负实数或者实部为负数的复数，为此可对式（34）进一步简化为

[A1q≈2πts-1eqtμs-1！+                     4-1sij=1N（1+xj）n-2Aje-1-xj1+xjt1+q+q-1xjn，] （35）

式中：[μ=Imq/Imq]。

由式（35）可知，本文所提方法将无穷积分转化为高斯-雅克比的积分，相对于文献[14]和文献[16]具有简洁性。

2.2 A2（q）、A3（q）和A4（q）的近似解

由式（8）及式（34），则A2（q）表示为

[A2q≈（-1）n-1ImiqReiq-Re-iq×2πtn-1e-qtn-1！-4ij=1N（1+xj）n-2Aje-1-xj1+xjt1-q-q+1xjn.] （36）

由于工程上的q为负实数或者实部为负数的复数，为此可将式（36）进一步简化为

[A2q=-4ij=1N（1+xj）n-2Aje-1-xj1+xjt1-q-q+1xjn]. （37）

由式（11）及式（37），则当q为负实数或者实部为负数的复数，则A3（q）表示为

[A3q=4ij=1N（1+xj）n-2Aje-1-xj1+xjt1-q-q+1xjn]. （38）

由式（15）及式（35），则当q为负实数或者实部为负数的复数，则A4（q）表示为

[A4q≈-2πtn-1eqtμn-1！+4-1n-1×                     ij=1N（1+xj）n-2Aje-1-xj1+xjt1+q+q-1xjn.] （39）

至此，获得了当q为负实数或者实部为负数的复数，则得到A1（q）—A4（q）的近似解。

3 算例

3.1 计算方法精度及效率验证

为验证本文所提方法的精度和分析效率，利用传统的数值积分进行了对比验证。数值积分计算A1（q）—A4（q）的计算表达式为

[A1q=12k=0ωu/Δω-1ei×k×Δωti×k×Δω-qn-e-i×k×Δωt-i×k×Δω-qn+ei×k+1×Δωti×k+1×Δω-qn-e-i×k+1×Δωt-i×k+1×Δω-qnΔω]， （40）

[A2q=12k=0ωu/Δω-1e-i×k×Δωti×k×Δω-qn-ei×k×Δωt-i×k×Δω-qn+e-i×k+1×Δωti×k+1×Δω-qn-ei×k+1×Δωt-i×k+1×Δω-qnΔω]， （41）

[A3q=12k=0ωu/Δω-1ei×k×Δωt-i×k×Δω-qn-e-i×k×Δωti×k×Δω-qn+ei×k+1×Δωt-i×k+1×Δω-qn-e-i×k+1×Δωti×k+1×Δω-qnΔω]， （42）

[A4q=12k=0ωu/Δω-1e-i×k×Δωt-i×k×Δω-qn-ei×k×Δωti×k×Δω-qn+e-i×k+1×Δωt-i×k+1×Δω-qn-ei×k+1×Δωti×k+1×Δω-qnΔω]， （43）

式中：[ωu、Δω]分别为积分上限和积分步长。

从式（40）—式（43）可知，采用数值方法需要给定积分上限和积分间距。为验证所提方法的正确性，A1（q）—A4（q）中的参数为n、q及t，故本算例的参数取值为：1）q=-50±100i，t=0.1 s， n=9和2） q=-300±400i， t=10.0 s， n=9。分别从积分上限和积分步长来研究本文方法的正确性，本文方法取30个高斯积分节点。

3.1.1 基于积分上限的验证

经试算积分步长[Δω]取10-4 rad/s，3种工况：工况1—工况3积分上限[ωu]取102、103、104 rad/s，计算结果如表1—表4。工况3的耗时为619.640 s，本文方法的耗时为0.032 s，本文方法具有较高的计算效率。

从表1—表4可知，当积分步长[Δω]=10-4 rad/s、q=-50±100i时，计算A1（q）—A4（q）时积分上限在3种工况下均一致，且与本文方法误差很小，说明[ωu]=102 rad/s时计算结果即可获得精确解；当q=-300±400i时，计算A1（q）—A4（q）时3种工况积分上限均不同，但当[ωu]=104 rad/s时最接近本文方法。

综上所述，本文方法在计算A1（q）—A4（q）时不受积分上限的影响，具有精度高的特点，也说明了本文方法的正确性。

3.1.2 基于积分步长的验证

通过3.1.1可知，积分上限[ωu]=104 rad/s时，计算q=-50±100i和q=-300±400i可获得高精度解。为此，本算例的计算基于积分步长的验证[ωu]取104 rad/s。3种工况积分步长[Δω]取：10-2、10-3、10-4 rad/s，计算结果如表5—表8所示。

由表5—表8可知，在积分上限[ωu]取104 rad/s、q=-50±100i时，计算A1（q）—A4（q）积分上限[Δω=10-3]和[10-4 rad/s]的计算结果一致且与本文方法比较接近，说明获得了精确解；在q=-300±400i时，计算A1（q）—A4（q）时3种工况积分上限均不同，但当[Δω]=10-4 rad/s时最接近本文方法。综上所述，本文方法在计算A1（q）—A4（q）时不受积步长的影响，具有精度高的特点，也说明了本文方法的正确性。

3.2 计算精度与高斯-雅克比积分节点数的关系

高斯系列积分的精度受积分节点个数的影响，积分节点个数越大，精度越高；为此本文研究了高斯积分节点个数对计算精度的影响。根据3.1节的结论，本文方法在计算参数分别为1） q=-50+100i、t=0.1 s， n=9和2）q=-300+400i、t=10.0 s、n=10时计算精度高，故下面的研究取上述参数进行。高斯积分节点数取5～30个。图2给出了q=-50+100i、t=0.1 s、n=9的积分节点数对计算结果的误差分析。图3给出了q=-300+400i、t=10.0 s、n=9的积分节点数对计算结果的误差分析。

由图2可知，在q=-50+100i、t=0.1 s、n=9时，在计算A2（q）和A3（q）时，积分节点数较少时精度较差，且实部的影响比虚部更大；而计算A1（q）和A4（q）时计算精度受高斯积分节点数影响不大。由图3可知，q=-300+400i、t=10.0 s、n=10时，计算所有的积分式时精度受高斯积分节点个数影响不大。对比图2和图3，高斯积分节点的个数受q的影响会有所不同，q的实部或虚部绝对值越小受高斯积分节点个数的影响越大；高斯积分节点的个数取25个时可具有较好的精度。

4 结论

本文针对振动领域里4个含有符号函数[sgn（ω）egiωtkiω-q-n]的无穷积分的表达式求解方法存在的问题，提出了一种高效算法，通过算例分析获得如下结论：

1）利用本文方法分析4个积分式的解时，无需考虑积分间距和积分上限的影响。而使用传统数值方法计算时，则受积分间距和积分上限的影响，且不同q值的积分解的精度受积分间距和积分上限的影响也会不同。因此，采用传统数值积分时可能会出现偏差很大的解，这一点需要引起注意。

2）本文所提解法的精度受高斯-雅克比积分节点个数的影响较大，q的实部或虚部绝对值越小受高斯积分节点个数的影响越大。从论文所获得结果来看，当取25个高斯积分节点，积分精度可不受q值实部或虚部的影响。为此，建议采用本文方法计算此类含有符号函数积分式的解时取25个高斯积分节点。

参考文献

[1] 衡楠，郑锋华，李金鹏．直升机机载扫描反射镜组件随机振动仿真分析[J]．机械设计与制造，2020（6）：290-293，297.

[2] 李创第，杨雪峰，李宇翔，等．组合体结构随机地震动响应的简明解法[J]．地震工程与工程振动，2023，43（1）：24-33.

[3] 葛新广，张梦丹，龚景海，等．频响函数二次正交法在Davenport风速谱下结构系列响应简明封闭解的应用研究[J]．振动与冲击，2021，40（21）：207-214.

[4] 谢娟娟，李晋，田震，等．考虑路面不平顺随机性的汽车过桥动力响应分析[J]．振动与冲击，2021，40（14）：299-306.

[5] 李创第，李宇翔，葛新广．路面不平顺激励下车辆系统振动响应谱矩的简明封闭解[J]．广西科技大学学报，2021，32（3）：13-20.

[6] 李创第，贺王涛，葛新广.卷积型非粘滞阻尼结构随机地震动系列响应求解的虚拟激励法[J].广西科技大学学报，2021，32（1）：78-84.

[7] 赵中伟，张永高.基于本征-虚拟激励法的大跨钢结构风振响应分析[J].空间结构，2020，26（1）：15-23.

[8] 林家浩，张亚辉.随机振动的虚拟激励法[M].北京：科学出版社，2004.

[9] 林家浩，张亚辉，赵岩.虚拟激励法在国内外工程界的应用回顾与展望[J].应用数学和力学，2017，38（1）：1-32.

[10] ALDERUCCI T，MUSCOLINO G．Fully nonstationary analysis of linear structural systems subjected to multicorrelated stochastic excitations[J]．ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems，Part A：Civil Engineering，2016，2（2）：C4015007.

[11] 张文首，林家浩.大跨度结构考虑行波效应时平稳随机地震响应的闭合解[J].固体力学学报，2004，25（4）：446-450.

[12] 何彬彬，冯磊，冯玉林，等.不同轨道不平顺谱激励下高铁列车走行性分析[J].铁道工程学报，2022，39（7）：36-43.

[13] 桂水荣，陈水生，万水.基于路面一致激励车桥耦合非平稳随机振动分析[J].振动.测试与诊断，2018，38（5）：908-915，1077.

[14] 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（6）：04014173.

[15] BARBATO  M，VASTA  M．Closed-form solutions for the time-variant spectral characteristics of non-stationary random processes[J].Probabilistic Engineering Mechanics，2010，25（1）：9-17.

[16] 葛新广.分数导数型阻尼耗能结构地震动响应的分析方法研究[D].上海：上海交通大学，2022.

[17] 张志超，张亚辉，林家浩.车桥耦合系统非平稳随机振动分析的虚拟激励-精细积分法[J].工程力学，2008，25（11）：197-204.

[18] DIETHELM K．An improvement of a nonclassical numerical method for the computation of fractional derivatives[J]．Journal of Vibration and Acoustics，2009，131（1）：014502.

Efficient numerical solutions for a class of infinite

integrals with signed functions

GE Xinguang1， LU Jiakang1， ZHANG Lirong1， LUO Zhen*2

（1. School of Civil and Architecture Engineering， Liuzhou Institute of Technology，

Liuzhou 545616， China; 2. School of Civil and Architecture Engineering， Guangxi University of

Science and Technology， Liuzhou 545006， China）

Abstract： To study the dynamic response of engineering structures under some random excitation based on frequency domain method， it is necessary to solve four infinite integrals with signed functions for which， however， there is computational efficiency and accuracy problem. Firstly， the calculation of four integral equations was converted into that of one integral equation according to the integral arithmetic. Secondly， an efficient and concise numerical solution for the above integral equations was presented based on the residue theorem and Gaussian-Jacobian integration. Finally， through numerical examples， the effect of the integration upper limit and integration interval on the traditional method and that of the number of Gaussian-Jacobian integration nodes on the method in this paper were studied. The results show that the calculation results of the traditional method are greatly affected by the integration interval and integration upper limit. However， the method in this paper can achieve high accuracy and computational efficiency by taking 25 Gaussian-Jacobian integration nodes.

Keywords： infinite integral; signed functions; Gaussian-Jacobian integral; efficient numerical solution

（责任编辑：罗小芬）

收稿日期：2023-10-11；修回日期：2023-10-25

基金项目：国家自然科学基金项目（51368005）；柳州工学院高层次人才项目（202201）资助

第一作者：葛新广，博士，副教授，研究方向：土木工程结构抗震、抗风研究

*通信作者：罗臻，硕士，高级工程师，研究方向：土木工程结构抗震、抗风研究，E-mail：20507257@qq.com