馬秋敏,楊 洪,2,于亦飛
(1.江蘇海洋大學(xué) 理學(xué)院,江蘇 連云港 222005;2.江蘇省海洋資源開(kāi)發(fā)研究院(連云港),江蘇 連云港 222005)
碳循環(huán)描述了碳在整個(gè)地球化學(xué)循環(huán)中的交換過(guò)程,是自然界中物質(zhì)和能量循環(huán)的重要組成部分。碳循環(huán)的一部分包括化合物的分解,允許碳被更新并以其他形式使用。這部分過(guò)程的一個(gè)關(guān)鍵組成部分是植物材料和木質(zhì)纖維的分解,而分解木質(zhì)纖維的關(guān)鍵因素是真菌。2021年美國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽A題要求用數(shù)學(xué)建模的方法研究真菌分解木質(zhì)纖維素纖維過(guò)程及其相互作用問(wèn)題。本文擬采用修正Frisch法、最小二乘法、微分方程等研究真菌分解木質(zhì)纖維素纖維的過(guò)程,以及不同種類(lèi)真菌間的相互作用、在環(huán)境影響下的真菌生存能力。
基于題中給出的假設(shè)條件,以及影響真菌的外部影響因素(濕度、菌種間的影響、pH值等)和內(nèi)部影響因素(真菌的內(nèi)在性狀,如生長(zhǎng)速度、耐濕性等),并且假設(shè)存在下列相關(guān)系數(shù)a,b,c,d,構(gòu)建真菌分解速率y與影響因素間的多元線性回歸模型[1]:
其中,x1為真菌的耐濕性,x2為真菌生長(zhǎng)率,x3表示菌種間的影響,x4為pH值,y為真菌分解速率。
由賽題資料可知,真菌分解速率y和生長(zhǎng)率x2呈線性關(guān)系,真菌分解速率y與耐濕性x1呈正比關(guān)系。進(jìn)一步研究真菌耐濕性x1、真菌生長(zhǎng)率x2對(duì)真菌分解速率y的影響,根據(jù)生長(zhǎng)率x2與分解速率y2的數(shù)據(jù),利用Eviews軟件作x2與y2的散點(diǎn)圖如圖1所示。由圖1可見(jiàn),x2與y2近似呈對(duì)數(shù)關(guān)系。
圖1 真菌生長(zhǎng)率x2與分解速率y2關(guān)系
將真菌生長(zhǎng)率x2與真菌分解速率y2進(jìn)行線性回歸,運(yùn)用最小二乘法得到估計(jì)結(jié)果如圖2。
圖2 x2與ln y2最小二乘估計(jì)結(jié)果
由圖2可知,
其中,括號(hào)內(nèi)的數(shù)字是t值。R2=0.94說(shuō)明總離差平方和中有94%被樣本回歸直線解釋?zhuān)挥?%未被解釋??梢?jiàn)回歸方程擬合結(jié)果較好,且t檢驗(yàn)和F檢驗(yàn)均符合要求,回歸方程顯著。
將(2)式化簡(jiǎn)可得:
同理,得到耐濕性x1與分解速率y1的最小二乘估計(jì)結(jié)果如圖3。
圖3 x1與ln y1最小二乘估計(jì)結(jié)果
由圖3可得
其中括號(hào)內(nèi)的數(shù)字是t值。R2=0.97說(shuō)明總離差平方和中有97%被樣本回歸直線解釋?zhuān)挥?%未被解釋。表明回歸方程擬合結(jié)果較好,且t檢驗(yàn)和F檢驗(yàn)均符合要求,回歸方程顯著。
將(4)式化簡(jiǎn)可得:
同時(shí),在處理真菌耐濕性與生長(zhǎng)率關(guān)系時(shí),可發(fā)現(xiàn)不同真菌耐濕性和生長(zhǎng)率大致呈線性關(guān)系(這里的線性關(guān)系可看作廣義的二元關(guān)系,只要兩個(gè)變量直接有關(guān)聯(lián)即可)。因此,在保持真菌分解速率不變的情況下將兩種因素進(jìn)行合并,其擬合方程為:
將不同因素影響下的y1、y2進(jìn)行加總,得出合并處理后的真菌分解速率模型為:
(1)競(jìng)爭(zhēng)關(guān)系標(biāo)準(zhǔn)模型如下:
其中:x5、x6分別為A、B菌種在t時(shí)刻的數(shù)量;rA、rB為A、B菌種的固定增長(zhǎng)率;NA、NB是環(huán)境資源容許的菌種最大數(shù)量;σ1表示對(duì)于供養(yǎng)A菌種的食物量而言,單位數(shù)量B菌種(相對(duì)NB)的消耗為單位數(shù)量A菌種(相對(duì)NA)消耗的σ1倍;σ2表示對(duì)于供養(yǎng)B菌種的食物量而言,單位數(shù)量A菌種(相對(duì)NA)的消耗為單位數(shù)量B菌種(相對(duì)NB)消耗的σ2倍[2]。
(2)模型求解
對(duì)于兩種有競(jìng)爭(zhēng)關(guān)系的真菌而言,當(dāng)為了生存爭(zhēng)奪有限的同一食物來(lái)源和生活空間時(shí),競(jìng)爭(zhēng)力較強(qiáng)的種群往往達(dá)到環(huán)境允許的最大數(shù)量。任取有競(jìng)爭(zhēng)關(guān)系的A菌種和B菌種,其中,A菌種是生長(zhǎng)緩慢的真菌菌株,更適應(yīng)在變化的環(huán)境下生存和生長(zhǎng),而B(niǎo)菌種生長(zhǎng)速度較快,但在相同的環(huán)境變化中不那么茁壯。因此,A菌種固定增長(zhǎng)率rA小于B菌種固定增長(zhǎng)率rB,而A菌種的競(jìng)爭(zhēng)能力高于B菌種,即σ1>σ2;在環(huán)境適宜時(shí),取rA=0.2,rB=0.5;σ1=2,σ2=0.5,獲得A、B菌種的數(shù)量關(guān)系與相軌線如圖4。
由圖4可見(jiàn),最后數(shù)值穩(wěn)定于x6=0,x5=10,即A菌種達(dá)到最大數(shù)量,B菌種滅絕。
圖4 A菌種與B菌種數(shù)量關(guān)系及相軌線
環(huán)境變化導(dǎo)致了A、B菌種的自然增長(zhǎng)率的快速變化。取rA=0.4,rB=0.3,σ1=2,σ2=0.5,獲取A、B菌種的數(shù)量關(guān)系與相軌線如圖5。
圖5 自然增長(zhǎng)率變化后的A、B菌種數(shù)量關(guān)系與相軌線
任取C菌種和D菌種,假設(shè)其是捕食關(guān)系,其中C菌種是捕食者,D菌種是被捕食者。C菌種主要以木質(zhì)纖維中的營(yíng)養(yǎng)物質(zhì)和D菌種為食,且當(dāng)捕食D菌種到一定程度則不再捕食;D菌種主要以木質(zhì)纖維中營(yíng)養(yǎng)物質(zhì)為食。
利用Volterra食餌-捕食者模型來(lái)解決以下問(wèn)題:
其中:t為時(shí)間,x(t)、y(t)為D、C菌種在t時(shí)刻的數(shù)量;r為D菌種獨(dú)立生存時(shí)的增長(zhǎng)率;d為C菌種獨(dú)自存在時(shí)的死亡率;a表示C菌種掠取D菌種的能力;b表示D菌種對(duì)C菌種的供養(yǎng)能力。
C菌種主要以木質(zhì)纖維中營(yíng)養(yǎng)物質(zhì)和D菌種為食,當(dāng)環(huán)境快速波動(dòng)時(shí),仍可以D菌種為食生存;而D菌種主要以木質(zhì)纖維中營(yíng)養(yǎng)物質(zhì)為食,當(dāng)環(huán)境快速波動(dòng)時(shí),則極易死亡。
取C菌種初始值為150,D菌種初始值為300,當(dāng)環(huán)境適宜時(shí),取r=2,d=1,a=b=0.01代入上述方程作菌種數(shù)量關(guān)系和相軌線如圖6。
圖6 捕食關(guān)系的C、D菌種數(shù)量關(guān)系及相軌線
由圖6可見(jiàn):C、D兩菌種相互制約,當(dāng)C菌種逐漸增多時(shí),D菌種逐漸減少;當(dāng)C菌種增加到一定數(shù)量時(shí),由于種群內(nèi)部競(jìng)爭(zhēng)激烈,又導(dǎo)致其數(shù)量減少,從而使得D菌種天敵減少,進(jìn)而導(dǎo)致D菌種數(shù)量增加;當(dāng)D菌種增加到一定數(shù)量時(shí),種群間競(jìng)爭(zhēng)加劇,從而又導(dǎo)致D菌種數(shù)量減少。如此循環(huán),相互制約發(fā)展。
當(dāng)環(huán)境快速波動(dòng)時(shí),取r=2,d=1,a=b=0.05,代入方程得圖7。從圖7可以看到,C、D兩菌種仍滿足上述關(guān)系,且在環(huán)境惡劣時(shí),其數(shù)量變化較為顯著。
圖7 D菌對(duì)C菌供養(yǎng)能力變化后C、D菌種的數(shù)量關(guān)系及相軌線
對(duì)于相互依存的兩種真菌而言,相互依存而共生是其普遍現(xiàn)象。任取互利共生關(guān)系的E菌種和F菌種進(jìn)行研究,其中,E菌種可單獨(dú)存在,F(xiàn)菌種為E菌種提供食物,同時(shí)不能離開(kāi)E菌種單獨(dú)存活,且有助于E菌種的生長(zhǎng)。
建立模型如下:
其中:x7、x8為E、F菌種在t時(shí)刻的數(shù)量;rE、rF為E、F菌種固定增長(zhǎng)率;NE、NF是環(huán)境資源容許的E、F菌種最大數(shù)量;σ3表示對(duì)于供養(yǎng)E菌種的食物量而言,單位數(shù)量F菌種(相對(duì)NF)的消耗為單位數(shù)量E菌種(相對(duì)NE)消耗的σ3倍;σ4表示對(duì)于供養(yǎng)F菌種的食物量而言,單位數(shù)量E菌種(相對(duì)NE)的消耗為單位數(shù)量F菌種(相對(duì)NF)消耗的σ4倍。
在溫度適宜時(shí),取rE=1.8,rF=1.5;σ3=0.1,σ4=3,作菌種數(shù)量關(guān)系與相軌線如圖8。
圖8 互利共生關(guān)系的E、F菌種數(shù)量關(guān)系與相軌線
從圖8可看出,當(dāng)t趨于無(wú)窮大時(shí),兩種真菌的最大數(shù)量均增大且相互依存,在無(wú)外界因素劇烈變化的條件下,兩種真菌數(shù)量趨于穩(wěn)定。
當(dāng)環(huán)境快速波動(dòng)時(shí),取rE=1.8,rF=1.5;σ3=0.3,σ4=3,得到E菌和F菌的數(shù)量關(guān)系及相軌線如圖9。
圖9 σ3變化后E、F菌種數(shù)量關(guān)系及相軌線
從圖9可看出,在環(huán)境快速波動(dòng)時(shí),隨著σ3的增加,兩種真菌間共生率增大,表明其促進(jìn)作用也增大,且兩種真菌數(shù)量隨時(shí)間增加達(dá)到平衡位置的值也增大。
不同種類(lèi)真菌的耐濕性和生活的最佳溫度都和真菌生長(zhǎng)率有著密切關(guān)系,也決定了真菌最適宜的生存環(huán)境條件。針對(duì)這兩種因素建模,分析其最適宜的生存環(huán)境條件,最終得到真菌生長(zhǎng)率與真菌所處環(huán)境的溫度和耐濕性的關(guān)系式。
由賽題所給條件可得生長(zhǎng)率x2、溫度x3與耐濕性x1的數(shù)據(jù),利用Eviews軟件作x2與x1的散點(diǎn)圖如圖10。由圖10可見(jiàn),x2與x1近似呈對(duì)數(shù)關(guān)系。
圖10 x1與x2散點(diǎn)圖
將真菌生長(zhǎng)率x2與真菌分解速率x1進(jìn)行線性回歸,運(yùn)用最小二乘法得到估計(jì)結(jié)果如圖11。
圖11 ln x21與x1的最小二乘估計(jì)結(jié)果
由圖11可得
R2=0.99,Rˉ2=0.99,F(xiàn)=593.54,DW=1.69其中,括號(hào)內(nèi)的數(shù)字為t值。R2=0.99說(shuō)明總離差平方和中有99%被樣本回歸直線解釋?zhuān)挥?%未被解釋??梢?jiàn)回歸方程擬合結(jié)果較好,且t檢驗(yàn)和F檢驗(yàn)均符合要求,回歸方程顯著[3]。
將(11)式化簡(jiǎn)可得:
同理可得圖12的散點(diǎn)圖。由圖12可見(jiàn),x2與x3近似呈線性關(guān)系。
圖12 x2與x3散點(diǎn)圖
擬合得x22=-2.194952+0.397 525x3
由t檢驗(yàn)可知,常數(shù)項(xiàng)不顯著,故舍去。再次擬合得:
按上述步驟檢驗(yàn),得知回歸方程顯著。
綜上可得,真菌生長(zhǎng)率與溫度及耐濕性的關(guān)系式為:
假設(shè)溫度為5~35℃,且x3的系數(shù)為正,由式(15)可知,溫度對(duì)生長(zhǎng)率始終起促進(jìn)作用。為指數(shù)函數(shù),且當(dāng)-0.500434+2.755192x1>0時(shí),即x1>0.181633時(shí),為單調(diào)遞增函數(shù),即真菌的生長(zhǎng)率隨著溫度與耐濕性的增大而增加。反之,若-0.500 434+2.755 192x1<0,則為單調(diào)遞減,且當(dāng)x1→∞時(shí),x21→0,即真菌生長(zhǎng)率只與溫度有關(guān)。
分析表明,干旱地區(qū)與半干旱地區(qū)真菌的耐濕性最弱,且干旱地區(qū)與半干旱地區(qū)晝夜溫差較大,因此真菌生長(zhǎng)率較低。喬木林和熱帶雨林地區(qū)真菌的耐濕性最強(qiáng),因此可視為只受溫度的影響,其生長(zhǎng)率與溫度成正比關(guān)系。溫帶地區(qū)真菌耐濕性較強(qiáng),且溫度適宜,生長(zhǎng)率明顯高于干旱地區(qū)、半干旱地區(qū)及喬木林和熱帶雨林地區(qū),適宜真菌生長(zhǎng)[4]。