邱舒婷 張志祥

(海軍工程大學(xué)計(jì)算機(jī)工程系 武漢 430033)

解常微分方程程序的蛻變測(cè)試方法*

邱舒婷 張志祥

(海軍工程大學(xué)計(jì)算機(jī)工程系 武漢 430033)

數(shù)值計(jì)算方法被廣泛運(yùn)用于實(shí)際生活中,解微分方程問題則是數(shù)值計(jì)算問題中的一個(gè)重要分支,但由于求解微分方程的程序很難直接找到一個(gè)測(cè)試Oracle來驗(yàn)證程序的有效性,進(jìn)而提出用蛻變測(cè)試方法檢測(cè)解微分方程程序的有效性。通過對(duì)解常微分方程的典型程序測(cè)試的案例分析,提出了一種求解常微分方程程序的蛻變測(cè)試方法,最后進(jìn)行一個(gè)實(shí)例驗(yàn)證該方法的有效性。

數(shù)值計(jì)算;微分方程;蛻變測(cè)試;Oracle問題

(Department of Computer Engineering, Naval University of Engineering, Wuhan 430033)

Class Number TP311

1 引言

數(shù)值計(jì)算程序?qū)ξ覀兊娜粘Ｉ罘浅Ｖ匾?不僅廣泛運(yùn)用于各種理論學(xué)科中,還在動(dòng)力工程和醫(yī)學(xué)實(shí)踐中的任務(wù)關(guān)鍵型和安全關(guān)鍵型應(yīng)用中有所體現(xiàn)。

由于數(shù)值計(jì)算中一些不可避免的截?cái)嗾`差,取整誤差以及計(jì)算進(jìn)程中的傳播誤差會(huì)影響到最終結(jié)果,并且找不到一種確切的解法來解決這個(gè)問題,因此在很多情況下,測(cè)試人員很難甚至是不能得出程序的預(yù)期輸出,進(jìn)而無法將程序的輸出結(jié)果與預(yù)期結(jié)果相比較,這就是軟件測(cè)試中常說的Oracle問題[2]。據(jù)此,人們提出三種測(cè)試方法來為數(shù)值計(jì)算程序提供一個(gè)可選的Oracle[1]：

1) 靜態(tài)檢查,進(jìn)行代碼走查檢驗(yàn)程序。

2) 動(dòng)態(tài)執(zhí)行,比較不同的數(shù)值計(jì)算程序生成的結(jié)果。

3) 利用不同的參數(shù)設(shè)置和輸出結(jié)果的一致性來運(yùn)行該軟件。

運(yùn)用代碼走查的方法檢驗(yàn)程序只適用于程序代碼量并不多的情況,但在復(fù)雜程序中此法效率很低。即使可以檢測(cè)到一些數(shù)值解和實(shí)驗(yàn)結(jié)果之間的差異,也很難去驗(yàn)證它們是由于程序bug還是在物理建模中的錯(cuò)誤引起的。并且在實(shí)際情況中,不同的環(huán)境設(shè)置可能會(huì)導(dǎo)致不同的結(jié)果。為解決這種Oracle問題,Chen等提出了蛻變測(cè)試技術(shù)[1]。

在運(yùn)用蛻變測(cè)試技術(shù)解決數(shù)值計(jì)算軟件測(cè)試中遇到的Oracle問題上,Tse T H等利用網(wǎng)格精度的不同對(duì)偏微分方程程序的執(zhí)行結(jié)果有不同的影響構(gòu)造出一條蛻變關(guān)系,測(cè)試了一個(gè)狄利克雷邊界條件下的橢圓型偏微分方程程序[3],實(shí)驗(yàn)表明網(wǎng)格精度越高,暴露出程序錯(cuò)誤的可能性越大,但在不同網(wǎng)格精度中產(chǎn)生的相對(duì)誤差值也十分近似,因此適用范圍較為局限。

姚奕等利用蛻變測(cè)試技術(shù)求解直角坐標(biāo)系上的圖形面積來對(duì)面向整數(shù)錯(cuò)誤檢測(cè)的蛻變測(cè)試方法進(jìn)行實(shí)例研究[4]。測(cè)試發(fā)現(xiàn)由于整型錯(cuò)誤產(chǎn)生的隱錯(cuò),實(shí)驗(yàn)結(jié)果顯示基于蛻變關(guān)系的整型錯(cuò)誤檢測(cè)方法可檢測(cè)出平時(shí)發(fā)現(xiàn)不了的隱式非預(yù)期輸出,有效地提高了檢測(cè)整型錯(cuò)誤的效率[11]。

由此可見,相較于一般的測(cè)試方法,蛻變測(cè)試技術(shù)在數(shù)值計(jì)算軟件測(cè)試中能夠更有效地檢驗(yàn)復(fù)雜數(shù)值計(jì)算程序的正確性。

在數(shù)值計(jì)算問題的眾多分支中,微分方程數(shù)值解法占有重要的地位。它以數(shù)值代數(shù)和逼近論等學(xué)科為基礎(chǔ),并反過來推動(dòng)這些相關(guān)學(xué)科的發(fā)展,在科學(xué)計(jì)算、工程技術(shù)等領(lǐng)域有著極其廣泛的運(yùn)用[5]。但微分方程的數(shù)值解法根據(jù)微分方程的不同類別各不相同,且很難直接得出其預(yù)期結(jié)果,而蛻變測(cè)試技術(shù)在測(cè)試微分方程數(shù)值計(jì)算程序中運(yùn)用還不是很廣泛,本文將研究用蛻變測(cè)試(Metamorphic Testing,MT)技術(shù)來緩解數(shù)值計(jì)算程序中的Oracle問題。

2 相關(guān)定義

2.1 蛻變測(cè)試相關(guān)定義

定義1 蛻變關(guān)系(Metamorphic Relation,MR)。假設(shè)程序P是用于計(jì)算函數(shù)f的,f：X→Y在定義域?yàn)閄,值域?yàn)閅上的函數(shù)。假設(shè)f滿足關(guān)系Rf,且關(guān)系Rf滿足當(dāng)f的n組輸入為x1,x2,…,xn∈X(n>1)時(shí),對(duì)應(yīng)的一系列函數(shù)輸入結(jié)果分別為f(x1),f(x2),…,f(xn)∈Y。這種關(guān)系Rf稱為關(guān)于f的一條蛻變關(guān)系。因此,若程序P是正確的,那么它一定滿足推導(dǎo)式(1)：

R(I1,I2,…,In) =Rf(P(I1),P(I2),…,P(In))

(1)

式中：I1,I2,…,In是程序P對(duì)應(yīng)于x1,x2,…,xn的實(shí)際輸入,P(I1),P(I2),…,P(In)則是對(duì)應(yīng)的輸出結(jié)果。

以正切函數(shù)為例。對(duì)于任意兩個(gè)輸入x1和x2,滿足x2=π+x1時(shí),一定有tanx2=tanx1。這種關(guān)系即可視為一條蛻變關(guān)系：

Rtan：x2=π+x1→tanx2=tanx1

(2)

因此在進(jìn)行測(cè)試時(shí),可通過驗(yàn)證式(1)、(2)的成立與否來驗(yàn)證程序P的輸出結(jié)果正確與否。

定義2 蛻變測(cè)試(Metamorphic Testing,MT)。蛻變測(cè)試是一種利用未暴露程序錯(cuò)誤的成功測(cè)試用例衍生出后續(xù)測(cè)試用例,并用衍生測(cè)試用例驗(yàn)證是否滿足程序的某些必要屬性來判斷程序正確性的技術(shù)。

測(cè)試人員通常會(huì)先在原始程序中植入一個(gè)變異,再利用蛻變關(guān)系生成衍生測(cè)試用例進(jìn)行程序測(cè)試。若程序未支持在該蛻變關(guān)系下應(yīng)滿足的屬性,則說明程序存在變異[8]。

2.2 常微分方程的數(shù)值計(jì)算

定義3 龍格庫塔法(Runge-Kutta)：龍格庫塔法[6]是一種在工程上應(yīng)用廣泛的高精度單步算法,用于數(shù)值求解微分方程。對(duì)于一階精度的拉格朗日中值定理有：

對(duì)于微分方程：y′=f(x,y)即y(i+1)=y(i)+h*k1,k1=f(xi+yi)。

當(dāng)用點(diǎn)xi處的斜率近似值k1與右端點(diǎn)xi+1處的斜率k2的算術(shù)平均值作為平均斜率k*的近似值,那么就會(huì)得到二階精度的改進(jìn)拉格朗日中值定理：

y(i+1)=y(i)+[h*(k1+k2)/2]

k1=f(xi+yi)

k2=f(x(i)+h,y(i)+h*k1)

依次類推,如果在區(qū)間[xi,xi+1]內(nèi)多預(yù)估幾個(gè)點(diǎn)上的斜率值k1,k2,…,km,(m>1)并用它們的加權(quán)平均數(shù)作為平均斜率k*的近似值,顯然能構(gòu)造出具有很高精度的高階計(jì)算公式。經(jīng)數(shù)學(xué)推導(dǎo)與求解,可以得出四階龍格庫塔公式,也就是在工程中應(yīng)用廣泛的經(jīng)典龍格庫塔算法：

y(i+1)=y(i)+[h*(k1+2*k2+2*k3+k4)/6]
k2=f(x(i)+h/2,y(i)+h*k1/2)
k3=f(x(i)+h/2,y(i)+h*k2/2)
k4=f(x(i)+h,y(i)+h*k3)

2.3 常微分方程求解示例

在實(shí)際工作中,解決一項(xiàng)現(xiàn)實(shí)問題需要成千上萬個(gè)微分方程組,人工求解微分方程的通解是很不現(xiàn)實(shí)的,所以微分方程的數(shù)值解顯得尤為重要。

在此以一個(gè)一階非線性常微分方程(3)為例,給出初始條件當(dāng)x0=1時(shí),y0=1且x∈[1,6]。

(3)

假設(shè)用“四階龍格庫塔法”解上述常微分方程問題。在給定邊界條件和步長h的情況下,程序可以計(jì)算出每隔h步長后的x和y的值直到迭代結(jié)束[7]。為說明蛻變測(cè)試的有效性,我們會(huì)在解常微分方程的程序中植入一個(gè)變異。圖1是原始程序的關(guān)鍵代碼。

圖1 關(guān)鍵代碼圖

假設(shè)通過將上述程序中的第7行代碼:

k2=(*pfunc)(x0+h/2,y0+k1*h/2)

替換為：

k2=(*pfunc)(x0+h/2,y0-k1*h/2)

來植入一個(gè)變異[10],并將設(shè)計(jì)的該變異程序版本記為Mutant1。將上述原始和植入錯(cuò)誤的程序編譯執(zhí)行后,得到的運(yùn)行結(jié)果趨勢(shì)如圖2所示。

圖2 非線性原始和變異程序運(yùn)行結(jié)果圖

從上述原始和植入變異過后的程序結(jié)果中可以看出,兩個(gè)程序的輸出值幾乎重合,說明兩個(gè)程序的運(yùn)行結(jié)果十分相似。隨著x值的逐步增大,y值逐漸遞減,且最終的y100值變異的運(yùn)行結(jié)果和原始結(jié)果相差甚微。因此,當(dāng)這種問題出現(xiàn)在實(shí)際應(yīng)用中時(shí),我們就很難發(fā)現(xiàn)錯(cuò)誤的存在。

3 解一階非線性常微分方程的蛻變測(cè)試方法

由于無法通過一般運(yùn)算求得一階非線性常微分方程的通解,進(jìn)而無法構(gòu)造出程序的預(yù)期輸出,待測(cè)程序缺乏測(cè)試Oracle,因此直接檢測(cè)出程序中是否存在錯(cuò)誤是有一定難度的。

通常在一般的一階常微分方程求解問題中,討論一階非線性常微分方程u′=f(t,u)的絕對(duì)穩(wěn)定性*絕對(duì)穩(wěn)定性：不管在理論還是應(yīng)用方面,單步法和多步法都必須是穩(wěn)定的。但這種穩(wěn)定有兩方面的限制：一是要求h∈(0,h0)充分小,而實(shí)際用的h∈(0,h0)都是固定的；二是只允許初值有誤差,往后各步計(jì)算都精確,而實(shí)際計(jì)算時(shí)每步都可能有舍入誤差,為了控制這種誤差的增長,需對(duì)多步法提出進(jìn)一步要求,即絕對(duì)穩(wěn)定性。非常困難,因而通?？紤]它在解u臨近的線性化方程,即簡化成討論一種簡單的線性常微分方程u′=a*u的絕對(duì)穩(wěn)定性,其中a為實(shí)數(shù)或復(fù)數(shù)。

由此引申,由于在用蛻變測(cè)試檢測(cè)程序的有效性時(shí),檢測(cè)一階非線性常微分方程u′=f(t,u)的有效性難度很大,同樣可以考慮檢測(cè)u臨近的線性化方程u′=a*u的有效性。為解決上述無法檢測(cè)出原始程序中的變異問題,提出基于檢測(cè)一階非線性微分方程有效性的蛻變測(cè)試一般方法,如圖3所示。

圖3 基于一階非線性微分方程的蛻變測(cè)試一般方法

蛻變檢測(cè)一階非線性常微分方程的一般方法步驟描述如下：

1) 對(duì)于一個(gè)一般的一階非線性常微分方程y′=f(x,y),必定存在一個(gè)解該方程的一般方法Runge-Kuta,此處以四階龍格庫塔法為例,將方法Runge-Kuta對(duì)應(yīng)的解一階非線性常微分方程程序記為P。

2) 實(shí)例化一階非線性常微分方程y′=f(x,y),簡化后的方程為一階線性常微分方程,記為y′=a*y(a為實(shí)數(shù)或復(fù)數(shù))。方程y′=a*y對(duì)應(yīng)的Runge-Kuta法記為Runge-Kuta1,對(duì)應(yīng)的實(shí)現(xiàn)Runge-Kuta1法解一階線性常微分方程的程序記為P’。

3) 對(duì)一階線性常微分方程y′=a*y構(gòu)造蛻變關(guān)系MRi(i=1,2,3…),用MRi驗(yàn)證實(shí)例化的程序P’,若程序P’存在錯(cuò)誤,則可推斷原始程序P中存在錯(cuò)誤。

4 一階非線性常微分方程蛻變測(cè)試實(shí)例

4.1 實(shí)例化一階非線性常微分方程

對(duì)于一個(gè)解一階非線性微分方程y′=f(x,y)的程序P,直接檢測(cè)出其中是否存在錯(cuò)誤是有一定難度的。若能在其對(duì)應(yīng)的解實(shí)例化后的一階線性微分方程y′=g(x,y)的程序P′中檢測(cè)出變異的存在,則能證明原始程序P的算法是有錯(cuò)誤的。

因此,將原始的非線性函數(shù)f(x,y)=-x*y2實(shí)例化為一個(gè)線性函數(shù)g(x,y)=u*y,作為程序的初始函數(shù)。令u=-3,此時(shí)實(shí)例化后的新微分方程記為一階線性微分方程(4)

(4)

同樣在程序中植入變異Mutant1,將上述實(shí)例化后的解一階線性微分方程的原始和植入錯(cuò)誤的程序編譯執(zhí)行后,得到的運(yùn)行結(jié)果趨勢(shì)如圖5所示。

圖4 線性原始和變異運(yùn)行結(jié)果圖

遺憾的是,從以上兩個(gè)程序運(yùn)行結(jié)果中依舊只能發(fā)現(xiàn)與圖2即最初的常微分方程程序運(yùn)行結(jié)果類似的規(guī)律,即隨著x值的逐步增大,y值依舊在逐漸遞減,且最終的y100值在變異后程序中的運(yùn)行結(jié)果與原始程序的運(yùn)行結(jié)果均近似為0(由于程序中數(shù)值精度有限,計(jì)算結(jié)果存在一定的誤差)。進(jìn)而說明,通過將微分方程實(shí)例化到一個(gè)相對(duì)簡單的一階線性常微分方程上,是不足以驗(yàn)證原始的解一階非線性微分方程程序是否存在錯(cuò)誤的。

4.2 解一階線性微分方程程序的蛻變測(cè)試方法

4.2.1 蛻變關(guān)系的構(gòu)造

通過對(duì)一階線性微分方程(4)構(gòu)造蛻變關(guān)系來驗(yàn)證一階線性常微分方程程序的正確性,進(jìn)而分析原始程序是否存在錯(cuò)誤。一階線性常微分方程(4)可表達(dá)為

(5)

解方程(5)的程序記為被測(cè)程序P′。根據(jù)一階線性微分方程的一般數(shù)值求解方法,可得出式(5)的通解為

y=e3*e-3x,c=e3

(6)

根據(jù)式(6)的性質(zhì),可以構(gòu)造出以下蛻變關(guān)系,如圖5所示。

圖5 y=e3*e-3x函數(shù)圖像

在上圖中可以看出隨著x值的逐漸增大,y值是單調(diào)遞減的,且y值始終大于0。根據(jù)這一條性質(zhì),可以構(gòu)造蛻變關(guān)系1,記為MR1。則MR1可描述為

MR1：?i,j∈且iy(xj)>0。

當(dāng)x取它的相反數(shù)-x時(shí),此時(shí)對(duì)應(yīng)的函數(shù)為y=e3*e3x,函數(shù)y與y’關(guān)于y軸對(duì)稱,且存在恒等式y(tǒng)(xk)*y(-xk)=e6。根據(jù)這一條性質(zhì),可以構(gòu)造出蛻變關(guān)系2,記為MR2。則MR2可描述為

MR2：?k∈,y(xk)*y(-xk)=e6。

以上兩條蛻變關(guān)系MR1,MR2就是測(cè)試執(zhí)行被測(cè)程序P’的所用到的蛻變關(guān)系。若程序執(zhí)行結(jié)果未滿足以上的關(guān)系,就有充分的理由斷定被測(cè)程序P’存在錯(cuò)誤。

4.2.2 測(cè)試結(jié)果

根據(jù)MR1和MR2在Windows環(huán)境下使用VC++編譯器執(zhí)行變異程序,選取程序的部分運(yùn)行結(jié)果如表1所示。其中最后一列的yx*y-x為由蛻變關(guān)系MR2構(gòu)造的衍生測(cè)試用例。

從表1中可以看出,對(duì)于任意逐漸增大的x值,y值始終在不斷遞減且y>0。因此蛻變關(guān)系MR1檢測(cè)不到程序的任何錯(cuò)誤。而在最后一列yx*y-x中,yx*y-x的值很明顯不等于e6,因此利用蛻變關(guān)系MR2,檢測(cè)出了解一階線性常微分方程的程序存在異常,即蛻變關(guān)系MR2的檢錯(cuò)能力大大高于蛻變關(guān)系MR1。這是由于MR2的輸入變換約束力更強(qiáng),較MR1更為復(fù)雜,與“優(yōu)先考慮輸入的蛻變關(guān)系表達(dá)式較復(fù)雜一方”的原則一致[9]。實(shí)驗(yàn)證明原始的解一階非線性常微分方程的程序存在異常,程序版本Mutant1存在缺陷。

表1 測(cè)試執(zhí)行結(jié)果

5 結(jié)語

因?yàn)槿鄙俜治龇椒?大多數(shù)微分方程通常通過數(shù)值方法解決。由于如前面章節(jié)所說的種種原因,對(duì)于這樣的數(shù)值方法可能不存在測(cè)試Oracle。而將微分方程簡化并聚焦到一階線性常微分方程問題上,會(huì)有助于解決一些問題,但是還不足以暴露出一些細(xì)微的錯(cuò)誤,比如本文中所描述的錯(cuò)誤。相反的,我們?cè)诤喕姆匠痰幕A(chǔ)上確定了蛻變關(guān)系并進(jìn)行了蛻變測(cè)試,然后設(shè)法從程序中找出問題所在。

本文以一階線性常微分方程為例,總結(jié)了蛻變測(cè)試是如何幫助減少數(shù)值計(jì)算軟件測(cè)試中的Oracle問題,并通過實(shí)驗(yàn)證明利用蛻變關(guān)系的構(gòu)造可以檢測(cè)出微分方程程序中的一些細(xì)微問題,在未來的更深入的微分方程的應(yīng)用中我們可以得到進(jìn)一步的鉆研。

[1] Chen T Y,Cheung S C,Yiu S M.Metamorphi-c testing：A new approach for generating next test cases[R].HKUST-CS98-01.Hong Kong,1998.

[2] Chen T Y,Huang D H,Tse T H,et al.Case studies on the selection of useful relations in metamorphic testing[C]// Symposium on Software Engineering&Knowledge Engineering.Madrid,Spain,2004：569-583.

[3] Chen T,Feng T H,Tse.Metamorphic testing of programs on partial differential equations：a case study[C]// Proceedings of the 26th Annual International Computer Software and Applications Conference. Baltimore,Marriott,USA,2002：327-333.

[4] 姚奕,黃松,稽孟雨.面向整數(shù)錯(cuò)誤檢測(cè)的蛻變測(cè)試方法研究[J].計(jì)算機(jī)工程與科學(xué),2012：52-56.

[5] 李榮華,劉博.微分方程數(shù)值解法[M].第四版.北京：高等教育出版社,2009：38-41.

[6] Khodabin M,Rostami M.Mean square numerical solution of stochastic differential equations by fourth order Runge-Kutta method and its application in the electric circuits with noise[J]. Advances in Difference Equations,2015：3-19.

[7] Chapra S C and Canale R P.Numerical Methods for Engineers：with Software and Programming Applications[M].New York：McGraw-Hill,2002.

[8] Offutt A J,Rothermel G,Zapf C.An experime-ntal evaluation of selective mutation[C]// Proceedings of the Fifteenth International Conference on Software Engineering.IEEE Computer Society,Washington DC,2003：100-107.

[9] Dong G W,Nie C H,Xu B W,et al.An effective iterative metamorphic testing algorithm based on program path analysis[C]//Proceedings of the 7th Annual International Conference on Quality Software.IEEE Computer Society,Washington DC,2007：292-297.

[10] WU Peng,SHI Xiao-Chun,TANG Jiang-Jun,et al.Metamorphic Testing and Special Case Testing:A Case Study[J]. Journal of Software,2005：1210-1214.

[11] Manolache L I,Kourie D G.Software Testing Using Model Programs[J].Software Practice and Experience,2001,31(13)：1211-1236.

版 權(quán) 聲 明

本刊已許可萬方數(shù)據(jù)庫、中國學(xué)術(shù)期刊(光盤版)電子雜志社在中國知網(wǎng)及其系列數(shù)據(jù)庫等產(chǎn)品中以數(shù)字化方式復(fù)制、匯編、發(fā)行、信息網(wǎng)絡(luò)傳播本刊全文。著作權(quán)使用費(fèi)與本刊稿酬一并支付。作者向本刊提交文章發(fā)表的行為即視為同意我編輯部上述聲明。

《艦船電子工程》編輯部

Metamorphic Testing of Ordinary Differential Equations Program

QIU Shuting ZHANG Zhixiang

Numerical methods are widely used in our real life, solving differential equations is an important branch of numerical problems. But it is difficult to find a test Oracle to verify the effectiveness of programs which solving differential equations. In this paper, metamorphic testing method is used to detect the effectiveness of the program which used to solve differential equations. Through analyzing the case of a typical program case, we propose the metamorphic test method of solving ordinary differential equations. At last, an instance is verified to illustrate the validity of this method.

numerical calculation, differential equations, metamorphic testing, Oracle problem

2016年6月9日,

2016年7月25日

邱舒婷,女,碩士研究生,研究方向:軟件質(zhì)量保證技術(shù)。張志祥,男,副教授,碩士生導(dǎo)師,研究方向:程序設(shè)計(jì)方法,軟件工程等。

TP311

10.3969/j.issn.1672-9730.2016.12.007