雙腔光反饋干涉激光系統(tǒng)中Lang-Kobayashi方程的六階龍格-庫塔算法

于芳星， 姬 波， CHENG Quanrun， 盧紅星， 柳宏川

(1.鄭州大學(xué) 信息工程學(xué)院，河南 鄭州 450001； 2.伍倫貢大學(xué) 電子計(jì)算機(jī)與通信工程學(xué)院，澳大利亞 新南威爾士州 伍倫貢市 2522)

0 引言

由光反饋效應(yīng)逐漸形成的自混合(SMI)干涉理論被廣泛應(yīng)用于位移計(jì)算、震動(dòng)檢測、形貌重構(gòu)、粒子運(yùn)動(dòng)測量等領(lǐng)域[1]。Fischer等[2]首次提出光反饋干涉儀中的雙腔結(jié)構(gòu)(雙腔OFI系統(tǒng))由激光二極管(LD)和兩個(gè)外部高反射金鏡組成。該模型在研究具有許多自由度的非線性系統(tǒng)中的動(dòng)力學(xué)時(shí)，具有十分理想的效果。Yu等[3]和Fan等[4]先后研究了多重光反饋以及多模激光器的自混合干涉理論，提出了含預(yù)反饋的激光自混合干涉型位移測量結(jié)構(gòu)，可以在弱反饋條件下得到類鋸齒波以進(jìn)行判向。近年來，研究人員已在許多應(yīng)用中使用了雙腔OFI結(jié)構(gòu)：Jiang等[5]和Geng等[6]提出了一種基于雙腔OFI的多普勒速度測量方法； Zhang等[7]和Chen等[8]使用雙腔OFI實(shí)現(xiàn)了二維振動(dòng)測量和多個(gè)目標(biāo)的運(yùn)動(dòng)檢測；Mezzapesa等[9]提出了一種使用雙腔OFI來實(shí)現(xiàn)納米級位移傳感的方法。

Lang等[10]首次提出一種描述單模半導(dǎo)體激光器的速率方程：Lang-Kobasyashi(L-K)方程，它可以完整地描述具有外部光反饋(EOF)的單模激光的動(dòng)態(tài)行為，成為研究光反饋?zhàn)曰旌细缮嫘?yīng)動(dòng)態(tài)問題的最經(jīng)典方法之一。L-K方程本質(zhì)是一個(gè)常微分方程組，具有變系數(shù)、非線性的特征，因此無法獲得解析解(符號解)，需要對其進(jìn)行數(shù)值求解。目前，已有的研究對雙腔OFI系統(tǒng)中的L-K方程主要通過歐拉法[11]、梯形公式法[12]、改進(jìn)歐拉法[13]、四階龍格-庫塔法[14]進(jìn)行求解。但這些方法存在著求解精度不夠高、計(jì)算速度不夠快的問題，導(dǎo)致生成的光電信號不能夠十分準(zhǔn)確地描述激光的動(dòng)態(tài)行為，降低了系統(tǒng)的測量精度。因此，為了進(jìn)一步提高光電信號雙腔OFI系統(tǒng)的動(dòng)態(tài)行為求解精度，本文提出一種求解L-K方程的六階龍格-庫塔算法，通過選取更多的區(qū)間點(diǎn)計(jì)算積分曲線的斜率平均值，使其更接近于真實(shí)值以提高精度，并將其應(yīng)用于雙腔OFI系統(tǒng)的移動(dòng)物體運(yùn)動(dòng)檢測仿真軟件中進(jìn)行仿真實(shí)驗(yàn)。

1 雙腔OFI系統(tǒng)的L-K方程

如圖1所示，雙腔OFI系統(tǒng)由激光二極管(LD)、光電二極管(PD)、聚焦透鏡(Lens)、目標(biāo)1(Target-1)、目標(biāo)2(Target-2)組成。其中Taerget-1為被測量目標(biāo)，其表面的反射率非常低，LD與Target-1之間的腔稱為被測量腔。Target-2表面具有較高的反射率可以用于提供預(yù)反饋，LD與Target-2之間的腔稱為控制腔。

圖1 雙腔OFI系統(tǒng)Figure 1 Dual-cavity OFI system

L-K方程包括3個(gè)聯(lián)立的延遲微分方程(DDE)，通過修改L-K微分方程，可以擴(kuò)展得到雙腔L-K方程：

E(t-τ2)·cos(ω0τ2+φ(t)-φ(t-τ2))；

(1)

(2)

(3)

雙腔L-K方程中參數(shù)的物理意義和數(shù)值見表1。

表1 L-K方程中參數(shù)的物理意義Table 1 Physical meaning of the parameters in the L-K equation

2 雙腔OFI系統(tǒng)L-K方程的龍格-庫塔算法

2.1 龍格-庫塔算法

在自然科學(xué)和社會(huì)科學(xué)的許多領(lǐng)域，經(jīng)常會(huì)遇到一階微分方程的初值問題：

(4)

在一系列離散的點(diǎn)x1,x2,x3，…，xn上求出未知函數(shù)y(x)之值y(x1),y(x2),y(x3),…,y(xn)的近似值y1,y2,y3，…，yn，即為所求的微分方程初值問題的數(shù)值解[15]。龍格-庫塔算法是目前最為常見的一種微分方程初值問題的數(shù)值解法，相比于歐拉法、梯形公式法、改進(jìn)歐拉法具有較高的精度和穩(wěn)定性，而且容易編程，所以被廣泛應(yīng)用于各個(gè)領(lǐng)域。

龍格-庫塔算法的基本思想是在(xn,xn+1)之間取一些積分曲線的若干點(diǎn)的切線斜率，再進(jìn)行一次(或多次)算數(shù)(或加權(quán))平均后產(chǎn)生新的斜率，隨后按這個(gè)斜率從(xn,yn)以直線代曲線向前推進(jìn)一步。根據(jù)差商代替導(dǎo)數(shù)的思想，由微分中值定理可知：

(5)

注意到一階微分方程y′(x)=f(x,y) 可得

y(xn+1)=y(xn)+hf(xn+θh,y(xn+θh))。

(6)

2.2 四階龍格-庫塔算法

當(dāng)在區(qū)間內(nèi)取4個(gè)點(diǎn)計(jì)算平均斜率時(shí)，即為四階龍格-庫塔算法。四階龍格-庫塔公式并不唯一，其中最經(jīng)典的為[16]

(7)

2.3 六階龍格-庫塔算法

六階龍格-庫塔算法需要計(jì)算7個(gè)k值，根據(jù)上述推導(dǎo)過程，可以求得一種六階龍格-庫塔公式[17]：

(8)

2.4 L-K方程的龍格-庫塔算法

根據(jù)六階龍格-庫塔公式，由式(1)可以得到L-K方程的龍格-庫塔公式：

(9)

根據(jù)式(9)，可以得到L-K方程六階龍格-庫塔算法的偽代碼如下所示。

算法L-K方程的六階龍格-庫塔算法。

輸入：E、φ、N關(guān)于t的微分方程和初始值；

輸出：E、φ、N值。

① 設(shè)置E、φ、N的初始值E0、φ0、N0；

② 設(shè)置步長h；

③ for (k=1∶SamplingNumber) do

④ 計(jì)算微分方程中目標(biāo)的相位變化值；

⑤ 根據(jù)六階龍格-庫塔公式計(jì)算k1、k2、k3、k4、k5、k6、k7；

⑥ then

⑦ 依次計(jì)算E、φ、N；

⑧ 更新E0、φ0、N0；

⑨ end for

用n表示采樣數(shù)的值，那么六階龍格-庫塔算法的時(shí)間復(fù)雜度為O(n)。

3 仿真軟件設(shè)計(jì)和實(shí)現(xiàn)

雙腔OFI系統(tǒng)仿真軟件通過MATLAB進(jìn)行開發(fā)，用于測量移動(dòng)物體在運(yùn)動(dòng)中的動(dòng)態(tài)行為，觀察仿真結(jié)果是否符合預(yù)期并預(yù)判測量精度，最終指導(dǎo)硬件平臺上的實(shí)際物理測試實(shí)驗(yàn)過程。

雙腔OFI系統(tǒng)仿真軟件進(jìn)行模塊化后可以分為7個(gè)主要模塊，其模塊分解圖見圖2。

圖2 模塊分解圖Figure 2 Module diagram

雙腔OFI系統(tǒng)仿真軟件中的接口包括常量、變量、計(jì)算反饋強(qiáng)度、計(jì)算相位變化量、定義微分方程、龍格-庫塔算法、繪制圖像7個(gè)接口，接口定義示例見表2。

表2 常量接口Table 2 Constant interface

雙腔OFI系統(tǒng)仿真軟件的總體仿真流程為：通過六階龍格-庫塔算法求解L-K方程，然后根據(jù)所得的E(腔內(nèi)電場強(qiáng)度)值繪制OFI信號關(guān)于時(shí)間t的圖像。

4 仿真實(shí)驗(yàn)與結(jié)果分析

本實(shí)驗(yàn)的主要工作是根據(jù)電場強(qiáng)度E的值體現(xiàn)L-K方程的求解精度，不同反饋強(qiáng)度下的實(shí)驗(yàn)為對照組。因此，實(shí)驗(yàn)的目標(biāo)是觀察反饋強(qiáng)度對干涉效應(yīng)的影響及不同反饋強(qiáng)度下的測量精度。k1、k2的值由實(shí)驗(yàn)輸入?yún)?shù)反饋系數(shù)c1、c2求解得出。

4.1 單腔實(shí)驗(yàn)結(jié)果

在單腔實(shí)驗(yàn)中，設(shè)置c1=0，分別設(shè)置c2=0.3、c2=2、c2=3進(jìn)行3組實(shí)驗(yàn)。作為對照，實(shí)驗(yàn)分別采用歐拉法、四階龍格-庫塔算法(RK4)、六階龍格-庫塔算法(RK6)進(jìn)行求解。實(shí)驗(yàn)對比結(jié)果見圖3、圖4，圖中橫坐標(biāo)t代表時(shí)間，縱坐標(biāo)OFI signal代表系統(tǒng)生成的光反饋信號強(qiáng)度。

圖3 單腔實(shí)驗(yàn)中歐拉法與RK6的對比圖Figure 3 Comparison of Euler method and RK6 in single cavity experiment

圖4 單腔實(shí)驗(yàn)中RK4與RK6的對比圖Figure 4 Comparison of RK4 and RK6 in single cavity experiment

表3和表4展示了單腔實(shí)驗(yàn)中不同算法求解的部分E值和精度改進(jìn)率。通過計(jì)算可得，RK6相比于歐拉法的平均精度改進(jìn)率約為22.95%，RK6相比于RK4的平均精度改進(jìn)率約為6.01%。

對比實(shí)驗(yàn)結(jié)果分析表明，在歐拉法與RK6的對比實(shí)驗(yàn)中，當(dāng)c2較小時(shí)，其信號幅度也比較小，歐拉法和RK6的插補(bǔ)效果幾乎相同；當(dāng)c2分別增大到2和3時(shí)，歐拉法會(huì)出現(xiàn)較明顯的偏差，RK6的效果明顯優(yōu)于歐拉法。在RK4與RK6的對比實(shí)驗(yàn)中，當(dāng)c2較小時(shí)，RK4和RK6的插補(bǔ)效果幾乎相同；當(dāng)c2增大到3時(shí)，RK4會(huì)出現(xiàn)較明顯的偏差，RK6的效果明顯優(yōu)于RK4。

表3 單腔實(shí)驗(yàn)中歐拉法與RK6的部分E值Table 3 Some E values of Euler method and RK6 in single cavity experiment

表4 單腔實(shí)驗(yàn)中RK4與RK6的部分E值Table 4 Some E values of RK4 and RK6 in single cavity experiment

4.2 雙腔實(shí)驗(yàn)結(jié)果

在雙腔實(shí)驗(yàn)中，設(shè)置3組c1、c2的值分別為：c1=0.1，c2=0.3；c1=1.5，c2=0.3；c1=4，c2=0.3。實(shí)驗(yàn)對比結(jié)果見圖5、圖6。

圖5 雙腔實(shí)驗(yàn)中歐拉法與RK6的對比圖Figure 5 Comparison of Euler method and RK6 in dual-cavity experiment

圖6 雙腔實(shí)驗(yàn)中RK4與RK6的對比圖Figure 6 Comparison of RK4 and RK6 in dual-cavity experiment

表5和表6展示了在雙腔實(shí)驗(yàn)中不同算法求解的部分E值和算法精度改進(jìn)率。通過計(jì)算可得，RK6相比于歐拉法的平均精度改進(jìn)率約為22.94%，RK6相比于RK4的平均精度改進(jìn)率約為6.01%。

表5 雙腔實(shí)驗(yàn)中歐拉法與RK6的部分E值Table 5 Some E values of Euler method and RK6 in dual-cavity experiment

表6 雙腔實(shí)驗(yàn)中RK4與RK6的部分E值Table 6 Some E values of RK4 and RK6 in dual-cavity experiment

對比實(shí)驗(yàn)結(jié)果分析表明：在雙腔實(shí)驗(yàn)中，其總體趨勢與單腔實(shí)驗(yàn)相似，即RK6的效果明顯優(yōu)于歐拉法和RK4；但是當(dāng)c2=0.3，同時(shí)c1逐漸增大時(shí)，信號波形浮動(dòng)范圍進(jìn)一步增大，會(huì)導(dǎo)致歐拉法、RK4和RK6的插補(bǔ)效果都變得惡劣。

5 結(jié)論

本文通過分析微分方程數(shù)值解法的原理，選取更多區(qū)間點(diǎn)計(jì)算積分曲線的斜率平均值，提出一種求解雙腔OFI系統(tǒng)中L-K方程的六階龍格-庫塔算法，并在移動(dòng)物體運(yùn)動(dòng)檢測仿真軟件中進(jìn)行了對比實(shí)驗(yàn)。結(jié)果表明：六階龍格-庫塔算法與歐拉法相比，求解精度平均提高了約22%，與四階龍格-庫塔算法相比，求解精度平均提高了約6%。因此，通過六階龍格-庫塔算法求解L-K方程可以提高其求解精度，從而提高光電信號雙腔OFI系統(tǒng)對于移動(dòng)物體的高靈敏度感測精度。