基于MATLAB的二元線性回歸在小孔流速實(shí)驗(yàn)中的應(yīng)用

黃育紅，馬　杰，張慕華，王恒通，楊宗立

(陜西師范大學(xué)，陜西 西安　710119)

?

基于MATLAB的二元線性回歸在小孔流速實(shí)驗(yàn)中的應(yīng)用

黃育紅，馬杰，張慕華，王恒通，楊宗立

(陜西師范大學(xué)，陜西 西安710119)

摘 要：利用MATLAB軟件對(duì)相關(guān)物理量進(jìn)行二元線性回歸擬合，求出了擬合曲線的待定參數(shù)和置信區(qū)間，擬合曲線與實(shí)驗(yàn)測試數(shù)據(jù)點(diǎn)符合較好。總結(jié)出小孔流速的經(jīng)驗(yàn)公式，并對(duì)減少實(shí)驗(yàn)誤差的途徑進(jìn)行了討論。

關(guān)鍵詞：小孔流速；二元線性回歸；MATLAB；回歸模型；殘差及置信區(qū)間

水是最常見、最簡單的流體之一。由于水的粘滯系數(shù)較小[1]，常溫(10 ℃)時(shí)僅為1.307 mPa·s且不易壓縮，因此常被視為理想流。

1986年，蔡一坤等人[2]通過實(shí)驗(yàn)的方法得到小孔流速與同高度自由落體速度一樣的結(jié)論；無獨(dú)有偶，1995年，金進(jìn)生等人利用伯努利方程和流體運(yùn)動(dòng)的連續(xù)性方程推導(dǎo)出了上面的結(jié)論，并得出了小孔的流量率和液位下降的時(shí)間，即對(duì)圓筒形容器來說，液面下降時(shí)間跟液面高度成冪函數(shù)關(guān)系。

前人對(duì)小孔流速的研究都是從某一方面出發(fā)，綜合考慮孔徑的大小和液面高度對(duì)流水時(shí)間影響的文章未見報(bào)道，且都沒有總結(jié)出流水時(shí)間t關(guān)于孔徑d、水高h(yuǎn)的經(jīng)驗(yàn)公式。此外，二元線性回歸是大學(xué)物理實(shí)驗(yàn)中必須熟練掌握的知識(shí)，而小孔流速實(shí)驗(yàn)可以很好地體現(xiàn)二元線性回歸在處理數(shù)據(jù)中的實(shí)際應(yīng)用。文章選取6個(gè)底部開有不同直徑小孔的圓筒容器，系統(tǒng)地研究了小孔出口處的流水時(shí)間與孔徑、水高的關(guān)系。利用MATLAB軟件對(duì)相關(guān)物理量進(jìn)行二元線性回歸擬合，擬合出了流水時(shí)間t與孔徑d、水高h(yuǎn)的經(jīng)驗(yàn)公式，求出了擬合曲線的待定參數(shù)和置信區(qū)間，擬合曲線與實(shí)驗(yàn)測試數(shù)據(jù)點(diǎn)符合較好。最后討論了減小實(shí)驗(yàn)誤差的方法。

1實(shí)驗(yàn)過程及數(shù)據(jù)的MATLAB處理

文章主要利用MATLAB的編程和圖形繪制功能對(duì)容器排水時(shí)間進(jìn)行系統(tǒng)研究。其中函數(shù)庫中的cftool在線性回歸擬合中獨(dú)具優(yōu)勢，可以方便地進(jìn)行曲線擬合，并得到擬合曲線。同時(shí)編寫MATLAB語句，可分析擬合結(jié)果的方差、剔除實(shí)驗(yàn)壞值，得出置信區(qū)間等。通過R-square和RMSE(方均根)等參數(shù)評(píng)定擬合結(jié)果的優(yōu)劣和擬合度的高低，R-square越接近1，RMSE越小，擬合度越高，擬合結(jié)果越可信[6-7]。

選取6個(gè)規(guī)格相同的柱狀圓筒，其底部分別開有大小不同的圓孔，利用游標(biāo)卡尺對(duì)圓孔直徑進(jìn)行測量，得到直徑的平均值及標(biāo)準(zhǔn)偏差。每次實(shí)驗(yàn)向6個(gè)圓筒里注入相同高度的水，6次實(shí)驗(yàn)向筒內(nèi)依次注入不同高度的水，利用秒表測量在不同孔徑、不同水面高度情況下的流水時(shí)間。利用MATLAB分別繪制流水時(shí)間和水高的圖像(孔徑一定)，以及流水時(shí)間和孔徑的圖像(水高一定)。利用MATLAB的cftool命令對(duì)相關(guān)物理量的變化關(guān)系進(jìn)行二元線性回歸擬合，并給出殘差及置信區(qū)間，得到流水時(shí)間t關(guān)于孔徑d、水高h(yuǎn)的經(jīng)驗(yàn)公式。

1.1孔徑d的測量

表1　圓孔直徑d的實(shí)驗(yàn)測量值及標(biāo)準(zhǔn)偏差

1.2不同孔徑、不同水高時(shí)流水時(shí)間t的測量

為測量不同孔徑、不同水高時(shí)圓筒對(duì)應(yīng)的流水時(shí)間，我們?cè)O(shè)計(jì)6組實(shí)驗(yàn)，每組實(shí)驗(yàn)均向6個(gè)底部開有不同尺寸的圓孔的圓筒內(nèi)同時(shí)注入等高度的水，6組實(shí)驗(yàn)向圓筒注入水的高度依次為10 cm，15 cm，20 cm，25 cm，30 cm，40 cm。用秒表測量水從圓筒內(nèi)流完所用的時(shí)間t。圓孔直徑d、水面高度h和對(duì)應(yīng)的流水時(shí)間t列于表2。

表2　不同孔徑、不同水高的圓筒對(duì)應(yīng)的流水時(shí)間的測量

1.3利用MATLAB軟件繪制變量變化關(guān)系圖

利用MATLAB軟件，分別繪出孔徑d一定時(shí)，流水時(shí)間t關(guān)于水高h(yuǎn)的圖像(h-t圖線)和水高h(yuǎn)一定時(shí)，流水時(shí)間t關(guān)于孔徑d的圖像(d-t圖線)。圖1為流水時(shí)間隨水面高度的變化圖，其中實(shí)心菱形為實(shí)驗(yàn)測量點(diǎn)，實(shí)線為MATLAB擬合曲線，紅色線條表示筒徑為7.59 mm時(shí)流水時(shí)間隨水面高度的變化情況，其它彩色線條分別代表筒徑為11.52 mm，14.94 mm，19.11 mm，19.53 mm和23.25 mm時(shí)流水時(shí)間隨水面高度的變化情況。圖2為流水時(shí)間隨圓孔直徑的變化圖，其中實(shí)心菱形為實(shí)驗(yàn)測量點(diǎn)，實(shí)線為MATLAB擬合曲線，紅色線條表示水面高度為10 cm時(shí)流水時(shí)間隨孔徑的變化情況，其它彩色線條分別代表水面高度為15 cm，20 cm，25 cm，30 cm和40 cm時(shí)流水時(shí)間隨孔徑的變化情況。從圖1和圖2中均可以看出，當(dāng)圓筒底部的孔徑d不變時(shí)，h越大，t越大，即水面越高，流水時(shí)間越長，當(dāng)水面高度h不變時(shí)，d越大，t越小，即孔徑越大，流水時(shí)間越短。

圖1　流水時(shí)間t隨水高h(yuǎn)的變化

圖2　流水時(shí)間t隨孔徑d的變化

1.4小孔水流時(shí)間隨水高和孔徑變化的函數(shù)模型

利用MATLAB的cftool命令對(duì)圓筒流水時(shí)間關(guān)于不同孔徑、不同水高的函數(shù)關(guān)系進(jìn)行擬合。分別采用MATLAB中常見的7種函數(shù)模型對(duì)小孔流速問題進(jìn)行擬合，通過R-square和RMSE(方均根)等參數(shù)評(píng)定擬合結(jié)果的優(yōu)劣和擬合度的高低，R-square越接近1，RMSE越小，擬合度越高，擬合結(jié)果越可信。其中Power、Exponential和Polynomial模型擬合函數(shù)得出的R-square和RMSE值較為理想。利用上述三種模型擬合流水時(shí)間t隨不同水面高度h變化的函數(shù)曲線，擬合評(píng)定參數(shù)R-square和RMSE值如圖3所示。其中綠線、藍(lán)線和紅線分別代表Power、Exponential和Polynomial模型下擬和函數(shù)對(duì)應(yīng)的R-square和RMSE值。

圖3　三種模型下擬合h-t函數(shù)對(duì)應(yīng)的評(píng)定參數(shù)

利用Power、Exponential和Polynomial模型擬合流水時(shí)間t隨不同圓孔直徑d變化的函數(shù)曲線，擬合評(píng)定參數(shù)R-square和RMSE值如圖4所示。

圖4　三種模型下擬合d-t函數(shù)對(duì)應(yīng)的評(píng)定參數(shù)

其中綠線、藍(lán)線和紅線分別代表Power、Exponential和Polynomial模型下擬和函數(shù)對(duì)應(yīng)的R-square和RMSE值。從圖3和圖4中均可以看出，Power模型擬合出的R-Square值最接近于1，且相應(yīng)的RMSE值最小，因此Power模型是擬合小孔流速問題最佳的選擇，水流時(shí)間隨圓孔直徑和高度的變化關(guān)系t=f(h,d)可寫為函數(shù)t=khαdβ的形式，其中k、α、β為待定參數(shù)。

1.5利用MATLAB擬合函數(shù)表達(dá)式

利用MATLAB函數(shù)庫中的cftool工具確定水流時(shí)間隨圓孔直徑和高度的變化關(guān)系t=khαdβ中的待定參數(shù)，擬合得出的函數(shù)關(guān)系列于表3。為了方便處理數(shù)據(jù)，對(duì)函數(shù)t=khαdβ進(jìn)行線性化，兩邊取對(duì)數(shù)得到lnt=lnk+αlnh+βlnd，其中l(wèi)nd，lnh的相關(guān)數(shù)值見表3，對(duì)表2中列出的圓筒水流時(shí)間t隨不同孔徑d、不同水面高度h的測量值取對(duì)數(shù)，得到lnt的數(shù)值列于表4。

表3　擬合水流時(shí)間與高度和孔徑的函數(shù)及相關(guān)參數(shù)

表4　不同高度、不同孔徑對(duì)應(yīng)流水時(shí)間的對(duì)數(shù)值

根據(jù)表3和表4的數(shù)據(jù)，編寫如下的MATLAB程序[4]，擬合出最佳直線，求出待定參數(shù)，擬合出的相關(guān)物理量及評(píng)定參數(shù)列于表5。

MATLAB程序：

n=36；m=2;

Ltime=[2.511.66 1.27 1.15 0.83 0.42 2.81 1.921.52 1.46 1.11 0.72 3.00 2.11 1.72

1.661.29 1.19 3.15 2.26 1.88 1.79 1.74 1.34 3.26 2.38 1.98 1.88 1.83 1.45

3.452.55 2.16 2.08 2.06 1.64];

Ldiameter=[2.032.44 2.70 2.95 2.97 3.15 2.032.44 2.70 2.95 2.97 3.15 2.032.442.70 2.95 2.97 3.15 2.03 2.442.70 2.95 2.97 3.15 2.03 2.442.70 2.952.97 3.15 2.03 2.44 2.70 2.952.97 3.15];

Llength=[2.30 2.30 2.30 2.30 2.30 2.30 2.71 2.71 2.71 2.71 2.71 2.71 3.00 3.00 3.003.00 3.00 3.00 3.22 3.22 3.22 3.22 3.22 3.22 3.40 3.40 3.40 3.40 3.40 3.403.69 3.69 3.69 3.69 3.69 3.69];

A=[ones(n,1),Ldiameter ',Llength'];

[b,bint,r,rint,s]=regress(Ltime',A);

ssquare=sum(r.^2)/(n-m-1);

b,bint,s,ssquare

rcoplot(r,rint)

程序中n為樣本量，m為自變量數(shù)目，Ltime為lnt，Ldiameter為lnd，Llength為lnh，b表示回歸系數(shù)，程序運(yùn)行后給出的的b1、b2和b3分別代表lnk、β和α的數(shù)值，bint表示回歸系數(shù)的置信區(qū)間，r表示殘差，rint表示殘差的置信區(qū)間，s用于檢驗(yàn)回歸模型的統(tǒng)計(jì)量，程序運(yùn)行后有四個(gè)數(shù)值，分別為相關(guān)系數(shù)R2、F、F(1,n-2)分布中大于F值的概率P和殘差的方差s2，R2和F的數(shù)值越大，P和殘差的方差越小，說明回歸方程越顯著。

表5　MATLAB擬合出的參數(shù)值及置信區(qū)間

殘差及置信區(qū)間如圖5所示，每條紅線長度表示的是置信區(qū)間，黑色小圓圈代表殘差點(diǎn)，殘差圖中未出現(xiàn)異常點(diǎn)。從圖中可以看出，本次實(shí)驗(yàn)數(shù)據(jù)較好，沒有壞值，擬合出的結(jié)果較為可信。從而不同圓孔直徑排水時(shí)間t關(guān)于水高h(yuǎn)及孔徑d的經(jīng)驗(yàn)公式為t=e3.831 7h0.747 1d-1.568 5。

圖5　殘差及置信區(qū)間圖

2實(shí)驗(yàn)結(jié)果的討論分析

利用計(jì)算機(jī)軟件MATLAB處理極其繁雜的二元線性回歸問題，可使過程簡單清晰，結(jié)論更加明確，使物理問題更加直觀和形象。冪函數(shù)Power模型可較好地?cái)M合出流水時(shí)間關(guān)于孔徑d、水面高度h的變化關(guān)系。

從擬合函數(shù)t=e3.831 7h0.747 1d-1.568 5可以看出，孔徑不變時(shí)，水面越高，流水時(shí)間越長；水面高度不變時(shí)，孔徑越大，流水時(shí)間越短。但由于水不是理想的流體，實(shí)驗(yàn)者在實(shí)驗(yàn)測量過程及數(shù)據(jù)采集時(shí)會(huì)不可避免地引入誤差，實(shí)驗(yàn)者應(yīng)盡可能地減小誤差以獲得更加準(zhǔn)確的結(jié)果。實(shí)驗(yàn)中應(yīng)盡量地準(zhǔn)確讀數(shù)，可以采用多次測量取均值的辦法減小孔徑測量的誤差，并利用不確定度評(píng)定測量結(jié)果的準(zhǔn)確度。用秒表測量流水時(shí)間時(shí)，實(shí)驗(yàn)者及時(shí)迅速地操作秒表可減少誤差。此外，從實(shí)驗(yàn)曲線上可以看出，選取柱形圓筒的孔徑應(yīng)盡量不大于23.25 mm，因?yàn)榭讖皆酱?，流水速率越快，時(shí)間的測量就越不準(zhǔn)確。

文章將理論模擬和實(shí)驗(yàn)相結(jié)合，充分體現(xiàn)了MATLAB軟件在數(shù)值處理方面的優(yōu)勢，利用軟件進(jìn)行圖像繪制和函數(shù)擬合，使得物理規(guī)律清晰，物理結(jié)論可視化效果加強(qiáng)。

參考文獻(xiàn)：

[1]徐建峰，李榮勝，周靜，等．用布里淵散射測量水的粘滯系數(shù)[J]．光學(xué)學(xué)報(bào)，2001,9(21):1112-1115.

[2]蔡一坤，張建恒．測小孔流速實(shí)驗(yàn)[J]．物理實(shí)驗(yàn)，1986,6(6):244-246.

[3]金進(jìn)生，周祖利．小孔流速和流量率分析[J]．杭州大學(xué)學(xué)報(bào)，1995(22):90-9.

[4]徐金明.MATLAB實(shí)用教程[M]．北京：清華大學(xué)出版社，2005:2-3.

[5]鈔曦旭，楊萬民，唐純青．MATLAB及其在大學(xué)物理中的應(yīng)用[M]．西安：陜西師范大學(xué)出版社，2006.

[6]黃薇.用Matlab比較雙縫干涉和雙縫衍射[J].大學(xué)物理實(shí)驗(yàn)，2015(1)：4-5.

[7]鄭君剛.大學(xué)物理實(shí)驗(yàn)數(shù)據(jù)處理的Matlab應(yīng)用[J].大學(xué)物理實(shí)驗(yàn)，2015(2)：12-13.

The Application of Binary Linear Regression in the Experiment of Keyhole Flow Velocity Based on MATLAB

HUANG Yu-hong,MA Jie,ZHANG Mu-hua,WANG Heng-tong,YANG Zong-li

(Shaanxi Normal University,Shanxi Xi’an 710119)

Abstract：The relationships between water velocity at the exit of the keyhole and diameter as well as height of the water in the vessel have been systematically investigated.Dualistic linear regression fittings are also performed referring to the related physical quantities.The corresponding undetermined parameters and the confidence interval of the fitted curves are determined.The fitting curves are in good agreement with the experimental data.Furthermore,the empirical formula of the keyhole flow velocity is obtained and the method of reducing experimental errors is discussed.

Key words：keyhole flow velocity；dualistic linear regression；MATLAB；regression model；residual error and confidence interval

中圖分類號(hào)：O 4-39

文獻(xiàn)標(biāo)志碼：A

DOI：10.14139/j.cnki.cn22-1228.2016.001.024

文章編號(hào)：1007-2934(2016)01-0094-04

基金項(xiàng)目：國家自然科學(xué)基金(No.51302162)；陜西省自然科學(xué)基金(2014JQ1016)；陜西師范大學(xué)中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金資助(No.GK261001150)

收稿日期：2015-07-31