一類含符號函數的無窮積分的高效數值解法

葛新廣 盧嘉康 張梨榮 羅臻

摘 要：基于頻域法研究幾類隨機激勵下工程結構的動力響應，需要求解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