儲(chǔ)油罐的變位識(shí)別與罐容表標(biāo)定

/儲(chǔ)油罐的變位識(shí)別與罐容表標(biāo)定薛申芳（邢臺(tái)學(xué)院數(shù)學(xué)系，邢臺(tái):054０0１）摘要該文就儲(chǔ)油罐的變位識(shí)別與罐容表標(biāo)定問題1，用四次多項(xiàng)式灌內(nèi)油體離散變率進(jìn)行回歸得到連續(xù)變化率,建立了連續(xù)型微分方程模型,通過該微分方程初值問題，得到傾斜灌內(nèi)油體體積隨高度的變化規(guī)律，然后把無變位與傾斜情況的灌內(nèi)油體隨測(cè)量高度的變化引起的灌內(nèi)油量差異，以及當(dāng)灌內(nèi)油體相同時(shí)，而油標(biāo)高度測(cè)量值差差異來判斷分析影響。該文就問題2的傾斜旋轉(zhuǎn)油罐，通過坐標(biāo)系旋轉(zhuǎn)方法,建立了計(jì)算灌內(nèi)油體與油位關(guān)系的精確積分模型；由于灌內(nèi)油體不是油位、傾斜角以及旋轉(zhuǎn)角的初等函數(shù),利用所給測(cè)量數(shù)據(jù)和所建立的數(shù)學(xué)模型進(jìn)行變位參數(shù)的確定時(shí)采用了最小二乘方法，且對(duì)部分被積函數(shù)的二元二次多項(xiàng)式展開近似處理，以及對(duì)傾斜角和旋轉(zhuǎn)角的離散化（把傾斜角和旋轉(zhuǎn)角的范圍化為一些離散小區(qū)間）搜索,計(jì)算結(jié)果顯示傾斜角為，旋轉(zhuǎn)角為;最后利用所給出的模型和所求出的傾斜角和旋轉(zhuǎn)角給出灌容標(biāo)定值。關(guān)鍵詞數(shù)學(xué)建模；罐容表標(biāo)定;多項(xiàng)式回歸；坐標(biāo)旋轉(zhuǎn)變換;積分1問題通常加油站都有若干個(gè)儲(chǔ)存燃油的地下儲(chǔ)油罐,并且一般都有與之配套的“油位計(jì)量管理系統(tǒng)”,采用流量計(jì)和油位計(jì)來測(cè)量進(jìn)/出油量與罐內(nèi)油位高度等數(shù)據(jù)，通過預(yù)先標(biāo)定的罐容表(即罐內(nèi)油位高度與儲(chǔ)油量的對(duì)應(yīng)關(guān)系）進(jìn)行實(shí)時(shí)計(jì)算，以得到罐內(nèi)油位高度和儲(chǔ)油量的變化情況。許多儲(chǔ)油罐在使用一段時(shí)間后，由于地基變形等原因,使罐體的位置會(huì)發(fā)生縱向傾斜和橫向偏轉(zhuǎn)等變化(以下稱為變位)，從而導(dǎo)致罐容表發(fā)生改變。按照有關(guān)規(guī)定,需要定期對(duì)罐容表進(jìn)行重新標(biāo)定?，F(xiàn)用數(shù)學(xué)建模方法研究解決儲(chǔ)油罐的變位識(shí)別與罐容表標(biāo)定的問題。(1)對(duì)小橢圓型儲(chǔ)油罐（兩端平頭的橢圓柱體)，分別對(duì)罐體無變位和傾斜角為＝4。10的縱向變位兩種情況做了實(shí)驗(yàn)（實(shí)驗(yàn)數(shù)據(jù)參看CMCM2０10A附件1),去建立數(shù)學(xué)模型研究罐體變位后對(duì)罐容表的影響，并給出罐體變位后油位高度間隔為1cm的罐容表標(biāo)定值。(2）對(duì)中間部分為柱面，兩端為球冠面的實(shí)際儲(chǔ)油罐,建立罐體變位后標(biāo)定罐容表的數(shù)學(xué)模型，即罐內(nèi)儲(chǔ)油量與油位高度及變位參數(shù)(縱向傾斜角度和橫向偏轉(zhuǎn)角度）之間的一般關(guān)系、利用罐體變位后在進(jìn)/出油過程中的實(shí)際檢測(cè)數(shù)據(jù)（參看CMＣM２０1０A附件2）確定變位參數(shù),并給出罐體變位后油位高度間隔為10ｃm的罐容表標(biāo)定值.2問題1２.1問題1分析?圖2出油變化率圖2出油變化率圖1進(jìn)油變化率從附表1中，可以得到油位以及油體的改變量,從而可以建立微分方程模型。圖1和圖2分別給出了無變位時(shí)進(jìn)油和出油情況的罐內(nèi)油體增加量（出油時(shí)罐內(nèi)油體增加量取負(fù)）與油位增量之比(油體變化率)隨油位變的離散圖形。?在圖1中有幾個(gè)不正常的點(diǎn)，可能由于測(cè)量誤差造成的；而從圖2中的點(diǎn)比較正常,所以在無變位情況建模時(shí)采用附表1中的出油數(shù)據(jù)。?對(duì)傾斜情況,通過對(duì)附件1中的數(shù)據(jù)分析，采用了附件１中進(jìn)油測(cè)量數(shù)據(jù)，且在計(jì)算變化率時(shí)把一組點(diǎn):累計(jì)進(jìn)油量為0.0017ｌ(該值由積分計(jì)算得到)時(shí)，油位高度為0.上述兩種情況,通過對(duì)灌內(nèi)油體的離散變化率進(jìn)行擬合，得到連續(xù)變化率，進(jìn)而可以求解連續(xù)變量的微分方程去得到灌內(nèi)油體體積隨高度的變化規(guī)律。２．２問題1變量說明記：罐內(nèi)油體的體積（）；:罐內(nèi)油位()；：油位增量(）;：罐內(nèi)油體的增加量()；：，即為油體關(guān)于的變化率的擬合函數(shù);其中,為1時(shí)表示無變位情況,為2時(shí)表示傾斜情況.3問題1數(shù)學(xué)模型及其及解?下面通過建立模型，得到無變位和傾斜兩種情況下罐內(nèi)油體的體積與油位高度的關(guān)系。３.1無變位情況 由2。1的分析以及２。2定義的變量知，現(xiàn)選用出油數(shù)據(jù)對(duì)隨著油位高度變化的數(shù)據(jù)，用四次多項(xiàng)式進(jìn)行回歸得[1]：（1）這里顯著性水平取為，其決定系數(shù)為，,，可知回歸方程非常顯著，回歸模型成立。通過解微分方程初值問題：（2）得到,其解為（３)上式即為無變位情況下罐內(nèi)油體的體積與油位高度的關(guān)系.３.2傾斜情況為了讓結(jié)果更符合實(shí)際，利用附件１中傾斜進(jìn)油變位數(shù)據(jù),對(duì)隨著油位高度變化的數(shù)據(jù)用四次多項(xiàng)式進(jìn)行回歸得：(4）顯著性水平也取為,決定系數(shù)為，,,回歸模型成立（經(jīng)驗(yàn)證比用出油數(shù)據(jù)擬合效果好）。通過解微分方程初值問題:（５）得到，其解為(6)無變位與傾斜情況的灌內(nèi)油體隨測(cè)量高度的變化以及灌內(nèi)油量差參看圖3.當(dāng)灌內(nèi)油體相同時(shí),而油標(biāo)高度測(cè)量值差差異的部分?jǐn)?shù)據(jù)參看表1.傾斜灌容標(biāo)定值參看附表1。圖3無變位與傾斜比較圖3無變位與傾斜比較油量差表1相同油體對(duì)應(yīng)的油位Ｖ(m3)０。00１70．５０1７1.00１71．50172。0０１72。5０173。00173.５017h1(ｍ)0。00150．２1９0０。359４0。4８430．6０390。７２350.８４790。９870h２(m）０0。18５40.327０0.４5８20。58220．７０.8144０。９333注：上表中為無變位油位，為傾斜油位.４問題2?下面建立中間部分為柱面，兩端為球冠面的實(shí)際儲(chǔ)油罐變位后標(biāo)定罐容表的數(shù)學(xué)模型.4.1問題2分析 油罐變位包括傾斜和旋轉(zhuǎn)兩種變位方式，這里先對(duì)傾斜進(jìn)行討論,然后再加上旋轉(zhuǎn).首先建立適當(dāng)?shù)淖鴺?biāo)系，對(duì)于傾斜變位,相當(dāng)于實(shí)行坐標(biāo)系旋轉(zhuǎn),在新的坐標(biāo)系下，計(jì)算罐體中的油體與坐標(biāo)高度的關(guān)系;對(duì)于旋轉(zhuǎn)變位,根據(jù)罐體的軸對(duì)稱性，便可以得到油位的坐標(biāo)高度與測(cè)量浮標(biāo)高度的關(guān)系，從而得到灌內(nèi)油體與測(cè)油浮標(biāo)高度的關(guān)系的數(shù)學(xué)模型。然后利用最小二乘法得到參數(shù)和。4．２問題2變量記號(hào)記：油罐縱向傾斜角(弧度），假設(shè)不超過10o×π/18０），且假定油罐右端上斜;:油罐旋轉(zhuǎn)角(弧度，假設(shè)不超過10ｏ×π/180)；：油罐柱面半徑（1．5）;:油罐柱體的長度（8）；：油罐油面鉛垂坐標(biāo)（)；：變位后罐內(nèi)油浮子的高度,即浸泡在罐油中的浮桿長度(ｍ);V：罐油體積()；：球冠半徑（1．６25）。5問題２模型?現(xiàn)在建立問題2的數(shù)學(xué)模型.5。1無變位油罐的坐標(biāo)系對(duì)無變位油罐而言,以柱面左側(cè)圓的圓心為坐標(biāo)中心，指向正上方,軸穿過柱面中心軸指向右方，建立右手坐標(biāo)系(參看圖4）。圖4圖4罐體無變位坐標(biāo)系記分別為平面與罐體柱面所交矩形的四個(gè)頂點(diǎn)（參看圖４）.則在坐標(biāo)系下，油罐柱面方程為：（7）左、右球冠面滿足的方程分別為：（8)和（９）柱面左、右側(cè)所在的平面方程分別為：(1０）和(11）5。2傾斜油罐的坐標(biāo)系假設(shè)油罐右端向上傾斜,傾斜角度為，這等價(jià)于將坐標(biāo)系繞軸順時(shí)針旋轉(zhuǎn)角度,假設(shè)得到的新坐標(biāo)系為，則這兩個(gè)坐標(biāo)系之關(guān)系為［2]：或(12）其中為旋轉(zhuǎn)矩陣（正交矩陣),且(13)這里旋轉(zhuǎn)矩陣中的取旋轉(zhuǎn)角度的負(fù)值（繞軸逆時(shí)針旋轉(zhuǎn)取正，順時(shí)旋轉(zhuǎn)針取負(fù)）.于是在坐標(biāo)系下，油罐柱面方程為：(１4)左、右球冠面滿足的方程分別為：（15）和（16)柱面左、右側(cè)所在的平面方程分別為：（17)和（18)記分別為點(diǎn)在坐標(biāo)系下的坐標(biāo),則(19)且在假定最大傾斜角度不會(huì)超過10o的情況下,易得:(2０)５.3傾斜油罐灌容模型 下面在坐標(biāo)系下討論傾斜油罐內(nèi)油體與坐標(biāo)高度和傾斜角的關(guān)系：（２１)利用微元法知，在方向取微元，用過軸上的點(diǎn)作垂直于軸的平面去截取油罐體，得到的截面面積與有關(guān),記之為，則微元體的體積為,則罐油體積可用積分（2２）得到。由于罐面是的分片函數(shù)，故對(duì)的不同區(qū)間范圍進(jìn)行分別積分,即分為三個(gè)區(qū)間范圍分別進(jìn)行積分［3]，即(２3）其中為當(dāng)在不同區(qū)間上用過軸上的點(diǎn)作垂直于軸的平面去截取油罐體,得到的截面面積.分別表示用過軸上的點(diǎn)作垂直于軸的平面去截取油罐體左、右球冠部分所得到的截面面積。它們的表達(dá)式分別為：(24）(25)(26）（2５）（26）其中將和、的表達(dá)式代入(2３)式便得傾斜油罐內(nèi)油體與坐標(biāo)高度和傾斜角的關(guān)系的數(shù)學(xué)模型。5。４傾斜旋轉(zhuǎn)油罐灌容模型下面考慮加上旋轉(zhuǎn)發(fā)生（旋轉(zhuǎn)角度為,參看圖９)，為罐體傾斜和旋轉(zhuǎn)后的斜，且記：油罐最低點(diǎn)到油面的高度(參看圖5）;：油罐只發(fā)生傾斜時(shí)，游標(biāo)測(cè)得的高度（參看圖５和圖6）；圖5油罐傾斜后坐標(biāo)圖5油罐傾斜后坐標(biāo)圖6油罐旋轉(zhuǎn)斷面圖圖6油罐旋轉(zhuǎn)斷面圖由幾何關(guān)系以及所建的坐標(biāo)系，得:由于假定油罐右端往上傾斜，取負(fù)值，從而得到與的關(guān)系為（2７)將(27）式代入(2３）式，即得到了罐內(nèi)油體積與、及的關(guān)系的一般模型：（2８)至此，就建立了傾斜旋轉(zhuǎn)油罐灌容的精確數(shù)學(xué)模型.6問題2變位參數(shù)確定及灌容標(biāo)定（1）首先討論、的確定。利用附件２中E列(顯示油高)和Ｆ列（顯示油量容積）的的數(shù)據(jù)，即由測(cè)量值及得：(29）為確定和,理論上講課采用最小二乘法,令誤差為，即通過求極值問題：（3０)s.t．(３１）去得出和.遇到問題是并不是的初等函數(shù)。故采用離散方法和ＭATLAB軟件去解決。方法是把的區(qū)間［－１0π/18０，0]以及的區(qū)間［0，10π/180］均分成一些離散個(gè)小區(qū)間（每個(gè)小區(qū)間長度適當(dāng)小），對(duì)每個(gè)去得到,去求（３2)由于、的積分表示不是初等函數(shù),故采用近似處理，即把、積分表達(dá)式中的被積函數(shù)對(duì)、分別在(０，-0.６)、(0，8．3２）進(jìn)行二元二次多項(xiàng)式展開（帶來的體積誤差）,然后再進(jìn)行積分求得、的近似表示式,由于其解析式繁雜，用軟件進(jìn)行處理，不再寫出相應(yīng)的展開表達(dá)式。計(jì)算結(jié)果為(參看程序附錄2Vｐrog2_1).根據(jù)所得角度,把部分測(cè)量值相應(yīng)的理論值進(jìn)行了比較，絕對(duì)誤差（參看圖7）.圖圖7結(jié)果比較（2)罐容標(biāo)定。利用上述求得的角度值對(duì)罐容進(jìn)行標(biāo)定，計(jì)算結(jié)果參看圖8、計(jì)算程序參看附錄2Vprog2＿2［4].圖圖8標(biāo)定結(jié)果圖7總結(jié)?問題１用四次多項(xiàng)式灌內(nèi)油體離散變率進(jìn)行回歸得到連續(xù)變化率,建立了連續(xù)型微分方程模型，通過該微分方程初值問題,得到傾斜灌內(nèi)油體體積隨高度的變化規(guī)律，也可以采用問題2的方法去建立精確地積分模型得到傾斜灌內(nèi)油體體積隨高度的變化規(guī)律.問題2的精確地積分模型中，變位油罐內(nèi)油體體積隨高度的變化規(guī)律用積分表示的，由于不是油位、傾斜角以及旋轉(zhuǎn)角的初等函數(shù),在變?yōu)閰?shù)確定是帶來困難，計(jì)算機(jī)編程時(shí)只能近似處理，二元被積函數(shù)在何處進(jìn)行多項(xiàng)式展開是導(dǎo)致誤差大小的關(guān)鍵問題；為了克服灌內(nèi)油體體積的積分表達(dá)式難以表示成油位、傾斜角以及旋轉(zhuǎn)角的初等函數(shù),采用了把傾斜角和旋轉(zhuǎn)角的范圍化為一些離散小區(qū)間,在區(qū)間各分點(diǎn)得到測(cè)量數(shù)據(jù)與模型數(shù)據(jù)誤差平方和,然后搜索到使得該平方和最小的傾斜角以及旋轉(zhuǎn)角，該方法一會(huì)有一定的誤差，但傾斜角和旋轉(zhuǎn)角的離散化越細(xì)密誤差將會(huì)越小。??參考文獻(xiàn)[1]茆詩松，程依明,濮曉龍編。概率論與數(shù)理統(tǒng)計(jì)教程。北京：高等教育出版社,200４.7［２］劉林編。航天器軌道理論。北京：國防工業(yè)出版社,2000。6[３］趙樹嫄編。微積分。北京:中國人民大學(xué)出版社，2007.3［4]王沫然編.ＭAＴLAB與科學(xué)計(jì)算(第二版)。北京：電子工業(yè)出版社，2０06.7附錄1（問題1程序）%%無變位出油模型(油標(biāo)底高：0。０286m，對(duì)應(yīng)體積:0．0017l立方米，當(dāng)高度1.2時(shí)體積為２。6５23，灌體：4.1101ｍ３）%lｏadｆuｊian1％h=ｘｃ(:，1)；％v＝xc(：,2）;%dｈ=dｉｆf(h）;％dv=ｄiff（ｖ）；％dvｈ=—ｄｖ。／dｈ；％%ｐlot(ｈ（1：enｄ—１)，dvh，＇k.'）％[b,binｔ,ｒ，rｉｎｔ，ｓarｔs］＝regｒｅｓs(ｄvh,［ｈ(１：enｄ—1）．^４，h(1：ｅｎｄ-1）．^3,ｈ（１：enｄ-1）。＾2,ｈ（1：end-1），onｅs（sizｅ（h(1:end-1）)）]）;%%holdon％f=poly2sｙm（b）;%％ezplot(f，０，1.2）;％vh1=dｓolve('Dvh=-４57３48２79７70５335/5６29４99５３4２13１２＊h^4＋8０29３／764＊h＾３－8278６97６83８2２23/35１８4３72088８32*h^2＋7849５/5629499５34２131２*h＋2463485091７８9０97／２２5１799８1368５２48＇,＇vｈ（０）＝0',’h')；％％figure（2)％ezplot(vｈ1,０，1．2);hoｌdｏｎ；％ｑ1＝suｂｓ(vh１,１．２）％%變位進(jìn)油模型(油標(biāo)底高：0。０２86m,對(duì)應(yīng)體積:0.00１7l立方米，當(dāng)高度1。2時(shí)體積為2．6５23,灌體:４．110１m3）%loaｄfujｉａｎ1%ｈ=［0;xxj（:，1)]；％v=［0．001７1;ｘxj(：，２）］；％dh＝diff(h）;%dv=diff（ｖ)；%dvｈ=dv。/dｈ；％％plｏt（h（1:end－1），dvｈ,＇k：．＇）％［b,binｔ,r，rint,sａｒts]=regrｅｓs（ｄvh,[h（１：ｅnd—１）。^４,h（1：enｄ-1）。^3,h（1：end—１）。^2，h(1：end—１),ｏnｅs（size（ｈ（1:enｄ－1))）]）;%％holdon％ｆ=poly２sym(b）;%％ezpｌot（ｆ，0,１.２）;％vh2=dｓｏlve（’Dvh＝－７４62５81８６964704９／285６＊ｈ^4+73６１3/1428＊ｈ^３—53607／142８*h^2+4７3４295１６15９697／3518４37208８832＊h+812４４7５275880339/4５496＇，＇vh（０)=0．０017１',’ｈ＇）；%%fｉguｒe(２)%ｑ２＝subｓ(vh２,1。2)％ezplot(ｖｈ2,0，1.２）%ezplot（vh2-ｖh1）%axｉs([0，1.2，－1,5］)symsh％0。00170。５０171．0０1７1．５０１７2.００1７2。50１73．00173.5０17ｖ0=ｉnｐuｔ(＇v0＝')h01=doubｌe（solve（vh1-v0）)h02=ｄｏubｌｅ（ｓｏｌve（vh2－v0）)表1傾斜灌融表(標(biāo)定值）ｈvhvｈvｈv00.０0170。３０0０0.９023０。600０２．0７570.９１00３．40８20．０１00０。０20４0。31０00。938９0。61002。１1750.９20０3。４4870.02000.040４0。32００0。97570。6２０02。15950．93003．4８８70。０3000．０6160.３3001.01２70.6３００２．2017０.９4003.52810。0400０．08390。34０01．04980．64002．2４41０.950０３.56670。05００0。1０720。３５０01。０8７2０。65002.28660。96003.604７０。060０0.13160．３600１。12470．66002。32930。９70０3。6４18０.0７000。15700。37001．1６240。67002。３７220.９８0０3.6７８０0。08000．183２0.38001．2０020。69002。45840。9９003。7１３3０。090０0．２１０20．3９001．2３82０．70002．５０17１．０0003.747６0．10000.2３8１0.4０0０1。2764０.71002.54511.01０03.78０80。11000.26670.41001.31４8０。72002.58861.0２003。81270。12０00．２9590.420０１.35330．7３002。6３221．03003.843４0.１300０。３2590.4３001．39２００。74００２.67５9１。0400３.８7２80。１4000。３5６40。4４00１.430８０.75002。7１9７1。０50０３.９0０７0．15000。３８750。450０1．4６980.７6002.76３51．０６0０３．92700.１6000.41920．46001.50９０0。77002．80731.0７0０3。９5170.17000。４5１30。47001.54830。７800２。８51１1．0８０03。97470。1８000．48390。４８0０1。5878０。7９00２.8９501．０9003.９9５８0.19０00.51６９０.490０1.６2７50．80０0２．９３８71。1０0０4．01４90.20０0０。５５0４０.5０0０1.667３０.81００2.９8241．１1004．03200。210０0.58420.51０01。707３0．820０3。０２601．12０0４。０469０.22000.６１８4０．5２０01。74750.83003。06951．130０4．0５950.2３000．65290。53０0１。7８79０.84００３。11２81.14004.06960．24０00．68７70.5400１。8２850.8５00３．1５5９１.１５0０４。07７２0.25000。7２２9０。55001。86９２0．8６003.19８8１.16004.０82１0.２6000。75830.56001．9１０10．87003.2414１.17004.08410。27０00.79390。57００１.９5１20。8８０03.283７1。１80０4．0８3２0．280００．８2９80.58０01。99250。8９00３。325６1．１9０0４。０7910.290０0.8６600.5９００2.034０0。90003.3６711。20004．07１8附錄2（問題2程序)%Vｐｒog2_1（求α、β的MAＴＬAB程序)a＝1。5；%柱面半徑1．5z0=８；R=1.６25；%球體半徑t=-7＊pi/１80：０.01：0；t1＝0:0.01:6＊ｐi／180;ｌoadａlfａhk=x(1:４０：ｅｎd,１）；v＝x（1：40：end，2)；％容積測(cè)量值foｒｉ＝1:length（t）ｆorｊ＝１:ｌenｇth(t1）c=ｔ（ｉ)；%ｓin（t)近似地用t代替;d=１-ｔ（ｉ)^２/2；％ｃos(ｔ)近似地用1—t^2／２代替；ＸA=-3/2*ｄ;ＸB=－3/2＊ｄ—8＊c;XＣ＝3／2＊d－8＊c;XD=—XA;foｒk＝1:leｎｇth(hk）h=XA—2*c+1/d*（1．5＋（hk(ｋ）-１．5)*（1-t1（j）^2/2）)；ｓyｍsＸZ；f=sqｒt（R^2－（d*Ｘ+c*Z）^2－（-ｃ＊X+d＊Z－R＋1)＾2）；ft=ｔａy（ｆ，X，Z，0,—0.6）;％f的泰勒展開式（2階）,經(jīng)試驗(yàn)在z=－０.6處結(jié)果好（)，球冠理論體積為４。0５79，近似計(jì)算結(jié)果為4.0577Zm=ｄ*（Ｒ—1）－sqrt(ｄ^2*R^2—2＊d^2＊R+d^2—X^2—２＊ｃ＊X＊Ｒ+2*c＊X+２*R-１)；SL＝int(2＊ｆt,Z，Zm,c／d*X);％VL＝iｎt(SＬ，X,XＡ，ｈ）;ｇ=ｓqrｔ(Ｒ^2-(d*X+c*Z）^２—(-c＊X+d*Z-z０＋R－1）^２);gt=tay(g,X,Z,0,z0+0．32）;ZＭ=d＊（1＋ｚ0－R）+sｑｒt（d＾2*(1＋2＊z0-2＊R+z0^２—２＊z0＊R+R^２）-Ｘ＾２＋2＊ｃ＊X*（R-１－z0)—2*（z0-R—z0*R)－z0＾2—1）;SR=ｉnt（2*ｇt,Ｚ,z０/ｄ+c／d*Ｘ,ZM);%VR＝inｔ（ＳR,X，XＢ,ｈ）;ｆ１=sqｒt(a^2-（d＊X+c＊Z)＾2+1e-3）；％由于計(jì)算誤差，避免根式下出現(xiàn)負(fù)值，加了一個(gè)小的正數(shù).S1＝inｔ(2＊f1，Z，c/d＊X，—ａ／c-ｄ/c*X）;%V1=iｎt（S１，Ｘ，XＡ，h)+VL;S2=int（2＊f1,Z,c/d＊X，z0/ｄ＋c／ｄ*X)；％Ｖ2=int（S１,X，XA，XB)+int（S2，X，ＸB,h)+ＶL+ＶR;％柱體體積為三部分分別為：2.98，53．5682,56.５487。Ｓ3=ｉnt（２*f1,Z,a/c-d／c＊X，z0/d+c/d＊X);％V３=int（S1，X,ＸＡ,XＢ)+int(S2,X,ＸB,XD)+iｎｔ（S3，X，XD,ｈ）+int（SＬ,X，XA,XD)＋VR；ifh〈XAV(k)=real（int(S1，X，XA,h）+int（SＬ,X，XＡ，h));ｅlseｉfｈ>ＸB＆h＜＝XDV（ｋ)=real（iｎt(S1，X,XA，XB）+inｔ（S2,Ｘ，XB，h）+int(SL,X，XA，ｈ）+ｉnｔ(SR，X,XB，h）);elseＶ(k）=rｅal（int（S1，X,XA，XB）+int（Ｓ2,X，XB，ＸD)+int(S3,X,ＸD,h）+ｉnt（SL，X，XA,XD）+iｎt（SR,X,ＸＢ,h))；ｅndeｎｄS(i,j）＝sum（(v’-ｄouble(V))．^2）；i,jendenｄｍ=ｍiｎ（min(S)）;[U1,U2］=fｉnd(S==m）；alfa=ｔ（U1）ｂeｉtａ=t１（U2）％Vprog2_2（罐容標(biāo)定計(jì)算)a＝1。5；％柱面半徑１．5z0=8；R=1。62５;%球體半徑t=—0.０01；t１=0.1７1；lｏａdaｌfａhｋ=０：0.1:3；c=ｔ;％sin(t）近似地用t代替;ｄ=1-t^2／2；％cos（t)近似地用1—t^2/2代替;ＸA=-3/2＊d；XB=—3/2＊d-8＊c;XC=3／2＊d-8＊c；XD=－XＡ;ｆｏｒk=1:lenｇth（hｋ)h=ＸA-２*c+１/d*（１.5+（hk（k)-1．5