農(nóng)作物施肥效果分析

1、農(nóng)作物施肥效果分析摘 要我們通過研究氮、磷、鉀三種肥料對土豆和生菜的作用，來建立施肥量與產(chǎn)量關(guān)系的模型。通過回歸分析的方法，將所給的數(shù)據(jù)進行MATLAB工具箱擬合，并利用殘差分析的方法，建立反映施肥量與產(chǎn)量關(guān)系的模型并檢驗分析,找到產(chǎn)量的最優(yōu)解以及氮、磷、鉀三種肥料的最優(yōu)配合比，在耕地面積一定的情況下研究土豆或生菜可以達到得最大收益值。由此我們建立的土豆產(chǎn)量模型為w1=-12.8361+0.1903n+0.0842p+0.0735k-0.0003n2-0.0002p2-0.0001k2生菜產(chǎn)量模型為w2=-0.4938+0.0756n+0.0234p+0.0067k-0.0002n2求解得到土

2、豆產(chǎn)量的最值，當n=317.1667，p=210.5000，k=367.5000時，得出w1max=39.71，氮磷鉀肥料的最優(yōu)配合比為1.5:1:1.74，土豆是喜鉀作物。我們可以得出生菜的最值，當n=224，p=685，k=372時，得w2=24.53，可以看出生菜是喜磷作物。在應(yīng)用方面，為了直觀的展示最大的利潤以及最優(yōu)配合比，設(shè)計了一個GUI人機交互界面，這樣可以清晰明了表示獲得的最大收益值。關(guān)鍵詞：回歸分析 MATLAB擬合 殘差分析 最優(yōu)配合比 GUI人機交互界面 18一 問題重述俗話說“民以食為天”，我們的生活與農(nóng)作物的供應(yīng)息息相關(guān)。近年來，隨著人口增多，耕地減少，所以化肥對農(nóng)作物

3、的生長、提高農(nóng)作物的產(chǎn)量具有重要的意義。農(nóng)作物除了吸收水分和空氣中二氧化碳以獲得碳、氫、氧等元素外，還必須從土壤再吸收氮、磷、鉀和其他礦質(zhì)養(yǎng)分，并在太陽能的幫助下合成有機物質(zhì)，以建造自己的有機機體，但土壤中的常量營養(yǎng)元素氮、磷、鉀和其他礦質(zhì)養(yǎng)分一般不能滿足作物生長的需求，需要施用含氮、磷、鉀的化肥來補充。在本問題中，某研究所通過研究氮、磷、鉀三種肥料對土豆和生菜的作用，來建立施肥量與產(chǎn)量關(guān)系的模型。實驗中將每種肥料的施用量分為10個水平，在考察其中一種肥料的施用量與產(chǎn)量的關(guān)系時，把另兩種肥料固定在第7個水平上，通過回歸分析的方法，將所給的數(shù)據(jù)進行MATLAB擬合，從而建立反映施肥量與產(chǎn)量關(guān)系的

4、模型,找到產(chǎn)量的最優(yōu)解以及氮、磷、鉀三種肥料的最優(yōu)配合比，在耕地面積一定的情況下研究土豆和生菜可以達到的最大收益值，并從實際情況出發(fā)，評價該模型的優(yōu)缺點。二 問題假設(shè)及符號假設(shè)2.1假設(shè)在不同的實驗條件下，保持水分、溫度、濕度、光照、土壤狀況等外界條件一致。2.2 假設(shè)施加的化肥完全進入土壤，沒有揮發(fā)作用。2.3 假設(shè)在模型一、模型二中氮、磷、鉀三種肥料對土豆和生菜的影響是相互獨立的。2.4 符號假設(shè)：氮、磷、鉀施用量分別為n、p、k；土豆的產(chǎn)量為w1，生菜的產(chǎn)量為w2；氮肥每公斤價格為x，磷肥每公斤價格為y，鉀肥每公斤價格為z；土豆每噸利潤為a，生菜每噸利潤為b；其他的固定成本為m；總收益為

5、s。三 模型的建立和求解3.1 模型一的建立與求解為了大致分析土豆的產(chǎn)量w1、生菜的產(chǎn)量w2與n、p、k的關(guān)系，首先利用表中的數(shù)據(jù)分別做出w1與n、w1與p、w1與k的散點圖并進行擬合、建立模型，如下圖所示。圖一 w1對n的散點圖從圖一可以發(fā)現(xiàn)，隨著n的增加，w1先增加在減小，由MATLAB擬合為二次函數(shù)。其數(shù)學(xué)模型為w1=0+1n+2n2+；求參數(shù)可得w1=-0.0003395n2+0.1971n+14.74；求解最值得，當n=290.2798時，w1=43.3471；與實際最值得對比發(fā)現(xiàn)，模型給出的最值小于實際給出的最值，這是因為，給出的模型是為了探究一般情況，故取該曲線的最值點，而實驗給

6、出的最值點不在該曲線上，因而舍棄了實驗給出的最值點。圖二 w1對p的散點圖同理給出w1對p的散點圖，在擬合過程中我們發(fā)現(xiàn)指數(shù)函數(shù)擬合效果較二次函數(shù)擬合效果好，所以給出指數(shù)函數(shù)數(shù)學(xué)模型w1=45.99e-0.0002318p-13.78e-0.008642p求解函數(shù)最值得，當p=286.9502時，w1=41.8763，最值的誤差來源同樣是因為上述原因。圖三 w1對k的散點圖給出w1對k的散點圖，并得出數(shù)學(xué)模型w1=-1.449×10-9k4+2.257×10-6k3-0.001191k2+0.2603k+18.68由圖像可得最值，k=651時，w1=46.22。下面我們研究

7、生菜的產(chǎn)量w2與n、p、k的關(guān)系，并做出w2與n、w2與p和w2與k的散點圖并進行擬合、建立模型。圖四 w2對n的散點圖給出w2對n的散點圖，并得出數(shù)學(xué)模型w2=-0.0002381n2+0.1013n+10.23求解最值得，當n=212.7257時，w2=21.0046；最值的誤差是因為擬合方程與實驗數(shù)據(jù)不完全吻合。圖五 w2對p的散點圖給出w2對p的散點圖，并得出數(shù)學(xué)模型w2=-5.453×10-5p2+0.0606p+6.876求解最值得，當p=555.6574時，w2=23.7124，最值誤差來源仍是上述原因。圖六 w2對k的散點圖給出w2對k的散點圖，并得出數(shù)學(xué)模型w2=0

8、.004657k+16.27顯然由圖六所示，該模型為線性模型，故無最值的存在。3.2 模型二的建立與求解模型一只是給出了產(chǎn)量與單一變量（氮、磷、鉀的三種施用量）之間的關(guān)系，我們在此基礎(chǔ)上進行改進，建立產(chǎn)量與多變量之間的數(shù)學(xué)模型。首先先建立土豆產(chǎn)量與三種肥料施用量關(guān)系的模型。綜合上面的分析，土豆的產(chǎn)量與氮肥施用量的關(guān)系是二次函數(shù)模型；土豆的產(chǎn)量與磷肥施用量的關(guān)系是四次函數(shù)模型；土豆的產(chǎn)量與鉀肥施用量的關(guān)系是指數(shù)函數(shù)模型，為建立三者統(tǒng)一的回歸模型，必須簡化計算，將土豆的產(chǎn)量與磷肥施用量的關(guān)系、土豆的產(chǎn)量與鉀肥施用量的關(guān)系均轉(zhuǎn)化為二次函數(shù)模型。結(jié)合土豆產(chǎn)量的模型一建立如下回歸模型w1=0+1n+2p

9、+3k+4n2+5p2+6k2+其中0+1n+2p+3k+4n2+5p2+6k2是當給出氮肥施用量n、磷肥施用量p、鉀肥施用量k時，土豆產(chǎn)量w1的平均值，其中參數(shù)0，1，2，3，4，5，6稱為回歸系數(shù)，應(yīng)大致服從均值為0的正態(tài)分布。我們利用MATLAB統(tǒng)計工具箱中的命令regress求解，使用格式為b,bint,r,rint,stats=regress(w1,X，0.05)其中輸入w1為上述回歸模型中w1的數(shù)據(jù)（n維向量，n=30），X為對應(yīng)于回歸系數(shù)的=（0，1，2，3，4，5，6）的數(shù)據(jù)矩陣1 n p k n2 p2 k2（n*7矩陣，其中第一列為全1向量），0.05為置信水平；輸出b為的

10、估計值，bint為b的置信區(qū)間，r為殘差向量，rint為r的置信區(qū)間，stats為回歸模型的檢驗統(tǒng)計量，有4個值，第1個是回歸方程的決定系數(shù)R2（R是相關(guān)系數(shù)），第2個是F統(tǒng)計量的值，第3個是與F統(tǒng)計量對應(yīng)的概率值p，第4個是剩余方差s2。得到上述回歸模型的回歸系數(shù)估計值及其置信區(qū)間（置信水平0.05）、檢驗統(tǒng)計量R2、F、p、s2的結(jié)果見表參數(shù)參數(shù)估計值參數(shù)置信區(qū)間0-12.8361-20.6921，-4.980210.19030.1597，0.220820.08420.0418，0.126530.07350.0512，0.09584-0.0003-0.0004，-0.00035-0.000

11、2-0.0003，-0.00006-0.0001 -0.0001，-0.0000R2=0.9190 F=43.4925 p=0 s2= 6.1094表1由此可以得出該回歸模型為w1=-12.8361+0.1903n+0.0842p+0.0735k-0.0003n2-0.0002p2-0.0001k2同理，可得生菜與三種肥料施用量關(guān)系的模型二w2=0+1n+2p+3k+4n2+5p2+模型求解可得各參數(shù)的值參數(shù)參數(shù)估計值參數(shù)置信區(qū)間0-0.4938-5.0236，4.036010.07560.0468，0.104520.02340.0176，0.029230.00670.0006，0.01284

12、-0.0002-0.0003，-0.00015-0.0000-0.0000，-0.0000R2=0.7813 F=22.3293 p=0 s2= 4.4389表2由此可以得出該回歸模型為w2=-0.4938+0.0756n+0.0234p+0.0067k-0.0002n2四 模型的檢驗和分析在模型一的建立與求解過程中已對該模型進行了檢驗與分析，現(xiàn)在我們著重檢驗分析模型二。圖七 土豆產(chǎn)量殘差圖首先我們分析土豆產(chǎn)量的殘差圖，從圖七中發(fā)現(xiàn)在土豆產(chǎn)量的30個數(shù)據(jù)中，僅有2個數(shù)據(jù)（圖中紅線所示）的置信區(qū)間不包含零點，對于較多的點來說，殘差的置信區(qū)間均包含零點，這說明回歸模型能夠良好的符合原始數(shù)據(jù)。圖八

13、生菜產(chǎn)量殘差圖我們從圖八中發(fā)現(xiàn)在生菜產(chǎn)量的30個數(shù)據(jù)中，僅有3個數(shù)據(jù)（圖中紅線所示）的置信區(qū)間不包含零點，對于較多的點來說，殘差的置信區(qū)間均包含零點，這說明回歸模型能夠良好的符合原始數(shù)據(jù)。現(xiàn)在我們對模型二的結(jié)果進行分析，由表1顯示，R2=0.9190指因變量w1的91.90%可由模型確定且p=0，因而土豆產(chǎn)量的模型從整體來看是可用的。由表2顯示R2=0.7813指因變量w2的78.13%可由模型確定且p=0，因而生菜產(chǎn)量的模型從整體來看也是可用的。表1的回歸系數(shù)給出了土豆產(chǎn)量模型中的0，1，2，3，4，5，6的估計值，觀察他們的置信區(qū)間發(fā)現(xiàn)，均不含零點，所以說明各項均對模型來說是顯著的。表1的

14、回歸系數(shù)給出了生菜產(chǎn)量模型中的0，1，2，3，4，5，6的估計值，觀察他們的置信區(qū)間發(fā)現(xiàn)，除0外其他置信區(qū)間內(nèi)均不含零點，所以說明0次項對模型不顯著。用以下的MATLAB程序求出土豆產(chǎn)量的最值banana=(x)12.8361-0.1903*x(1)-0.0842*x(2)-0.0735*x(3)+0.0003*(x(1)2)+0.0002*(x(2)2)+0.0001*(x(3)2);x,fval=fminsearch(banana,0,0,0)求解得到最值，當n=317.1667，p=210.5000，k=367.5000時，得出w1max=39.71，由此結(jié)果我們可以得到氮磷鉀肥料的最優(yōu)

15、配合比為1.5:1:1.74。我們可以發(fā)現(xiàn)，w1max的值略小于實驗數(shù)據(jù)中給出的最大值，下面我們分析產(chǎn)生誤差的幾點原因。第一， 給出的實驗數(shù)據(jù)較少，不能很好的給出擬合模型。第二， 在模型二中，只考慮了氮、磷、鉀三種肥料獨立影響的情況，而忽略了三者之間的交互影響。第三， 在模型二中，缺乏正交試驗數(shù)據(jù)，無法給出較準確地模型。綜上所述，在模型二中w1max的值應(yīng)當較實際值偏小，但是我們可以看出土豆對鉀肥的喜好>氮肥的喜好>磷肥的喜好，土豆是喜鉀作物。在生菜的模型二中，產(chǎn)量與磷肥鉀肥成線性關(guān)系，因而無法求出最值，但我們可以看到產(chǎn)量與氮肥成二次關(guān)系，因而氮肥影響更加顯著，所以隨著氮肥的增加生

16、菜產(chǎn)量先增加后減少，但磷肥鉀肥卻成線性關(guān)系，雖然符合題中所給出的實驗數(shù)據(jù)，但不符合實際情況。因此我們只能由實際數(shù)據(jù)得到最值，當n=224，p=685，k=372時，得w2=24.53，可以看出生菜是喜磷作物。下面我們探究生菜產(chǎn)量的模型二與實際情況誤差較大的原因第一， 在給出磷肥鉀肥的實驗數(shù)據(jù)中，我們可以發(fā)現(xiàn)，產(chǎn)量基本是隨著施用量的增加而增加，因而擬合圖像只能是線性關(guān)系，對于這種情況，應(yīng)該增添磷肥鉀肥的實驗數(shù)據(jù)來完善模型二。第二， 在模型二中，只考慮了氮、磷、鉀三種肥料獨立影響的情況，而忽略了三者之間的交互影響。我們通過上述分析發(fā)現(xiàn)，土豆與生菜的模型二中，均未對交互項的影響做出探究，而氮磷鉀三者

17、的交互作用不可忽略，因此我們針對這一點對模型二進行改進。五 模型的改進在土豆產(chǎn)量的模型二中w1=0+1n+2p+3k+4n2+5p2+6k2+增加四項后得w1=0+1n+2p+3k+4n2+5p2+6k2+7np+8nk+9pk+10npk+我們把上述模型命名為模型三，在模型三中氮磷鉀三種肥料不是相互獨立的，其中，np為氮肥、磷肥之間的交互影響項，nk為氮肥、鉀肥之間的交互影響項，pk為磷肥、鉀肥之間的交互影響項，npk為氮肥、磷肥、鉀肥之間的交互影響項。模型三中的7、8、9、10稱為交互項影響系數(shù)，四個系數(shù)中的正值代表肥料的效果是相互促進的，下面直接給出生菜產(chǎn)量的模型三w2=0+1n+2p+

18、3k+4n2+5p2+6np+7nk+8pk+9npk+由于所給數(shù)據(jù)都是固定兩個水平，研究另一個的影響因素，題中沒有給出正交實驗數(shù)據(jù)，因而模型三中的交互項系數(shù)無法求出。但是模型三的優(yōu)點比較明顯，可以較準確地擬合出實際情況，并可以給出較準確的最值，以及可以比較幾個交互影響的權(quán)重，若想得出模型，建議完善實驗數(shù)據(jù)并給出多組正交實驗數(shù)據(jù)。六 模型的應(yīng)用價值現(xiàn)在我們主要分析模型二的應(yīng)用價值以及合理化建議，首先我們研究該模型的經(jīng)濟價值。土豆的經(jīng)濟收益模型是s=w1a-nx-py-kz-m生菜的經(jīng)濟收益模型是s=w2b-nx-py-kz-m我們可以根據(jù)產(chǎn)量和三種肥料之間的關(guān)系，某時間內(nèi)土豆生菜每噸的利潤，每

19、公斤化肥的成本以及固定成本等，利用上述兩個經(jīng)濟收益模型，計算出在一定耕地面積下，我們可以獲得的最大的土豆生菜種植的經(jīng)濟收益。為此我們設(shè)計了一個GUI人機交互界面來預(yù)測最大收入通過該模型可以指導(dǎo)我們合理使用化肥，降低對土壤酸度的影響，避免土壤板結(jié)、肥力下降，減少對土壤的污染，從而保護生態(tài)環(huán)境。七 參考文獻1姜啟源,謝金星，葉俊.數(shù)學(xué)模型.4版.高等教育出版社，2011.2陳桂明，戚紅雨，潘偉.MARLAB數(shù)理統(tǒng)計（6.x）.科學(xué)出版社，2002.3辛磊.如何做殘差分析. 八 附錄8.1 氮肥磷肥鉀肥作用下產(chǎn)量與施肥量散點圖的程序x1=0 34 67 101 135 202 259 336 404

20、 471;y1=15.18 21.36 25.72 32.29 34.03 39.45 43.15 43.46 40.83 30.75;x2=0 24 49 73 98 147 196 245 294 342;y2=33.46 32.47 36.06 37.96 41.04 40.09 41.26 42.17 40.36 42.73;x3=0 47 93 140 186 279 372 465 558 651;y3=18.98 27.35 34.86 38.52 38.44 37.73 38.43 43.87 42.77 46.22;x4=0 28 56 84 112 168 224 280

21、336 392;y4=11.02 12.70 14.56 16.27 17.75 22.59 21.63 19.34 16.12 14.11;x5=0 49 98 147 196 294 391 489 587 685;y5=6.39 9.48 12.46 14.33 17.10 21.94 22.64 21.34 22.07 24.53;x6=0 47 93 140 186 279 372 465 558 651;y6=15.75 16.76 16.89 16.24 17.56 19.20 17.97 15.84 20.11 19.40;%輸入六組數(shù)據(jù)scatter(x1,y1,'.

22、k')%繪出土豆在氮肥作用下產(chǎn)量與施肥量的散點圖grid onhold onp1=polyfit(x1,y1,2);%多項式擬合X1=1:0.001:1000;z1=polyval(p1,X1);plot(X1,z1)%輸出擬合圖像scatter(x2,y2,'.k')%繪出土豆在磷肥作用下產(chǎn)量與施肥量的散點圖grid onhold onp2=polyfit(x2,y2,2);%多項式擬合X2=1:0.001:1000;z2=polyval(p2,X2);plot(X2,z2)%輸出擬合圖像以下是應(yīng)用擬合工具箱指數(shù)擬合圖像scatter(x3,y3,'.k

23、9;)%繪出土豆在鉀肥作用下產(chǎn)量與施肥量的散點圖grid onhold onp3=polyfit(x3,y3,6);%多項式擬合X3=1:0.001:1000;z3=polyval(p3,X3);plot(X3,z3)%輸出擬合圖像scatter(x4,y4,'.k')%繪出生菜在氮肥作用下產(chǎn)量與施肥量的散點圖grid onhold onp4=polyfit(x4,y4,2);%多項式擬合X4=1:0.001:1000;z4=polyval(p4,X4);plot(X4,z4)%輸出擬合圖像scatter(x5,y5,'.k')%繪出生菜在磷肥作用下產(chǎn)量與施肥量

24、的散點圖grid onhold onp5=polyfit(x5,y5,5);%多項式擬合X5=1:0.001:1000;z5=polyval(p5,X5);plot(X5,z5)%輸出擬合圖像scatter(x6,y6,'.k')%繪出生菜在鉀肥作用下產(chǎn)量與施肥量的散點圖grid onhold onp6=polyfit(x6,y6,1);%多項式擬合X6=1:0.001:1000;z6=polyval(p6,X6);plot(X6,z6)%輸出擬合圖像8.2 三元函數(shù)二次擬合多項式及殘差圖程序w1=y1,y2,y3;n=x1,259*ones(1,20);p=196*ones(

25、1,10),x2,196*ones(1,10);k=372*ones(1,20),x3;X=ones(1,30);n;p;k;n.2;p.2;k.2;X=X'w1=w1'b,bint,r,rint,stats=regress(w1,X,0.05)%求出土豆產(chǎn)量與氮磷鉀施肥量三者函數(shù)式rcoplot(r,rint)%繪制殘差圖w2=y4,y5,y6;n=x4,224*ones(1,20);p=391*ones(1,10),x5,391*ones(1,10);k=372*ones(1,20),x6;X=ones(1,30);n;p;k;n.2;X=X'w2=w2'b

26、,bint,r,rint,stats=regress(w2,X,0.05)%求出生菜產(chǎn)量與氮磷鉀施肥量三者函數(shù)式rcoplot(r,rint)%繪制殘差圖8.3土豆生菜的N、P、K施肥總表格表1 土豆N、P、K施肥總表格產(chǎn)量N施肥量P施肥量K施肥量15363419637225.726719637232.2910119637234.0313519637239.4520219637243.1525919637213.4633619637240.8340419637230.7547119637233.46259037232.472592437236.062594937237

27、.962597337241.042599837240.0925914737241.2625919637242.1725924537240.3625929437242.7325934237218.98259196027.352591964734.862591969338.5225919614038.4425919618637.7325919627938.4325919637243.8725919646542.7725919655846.22259196651表2 生菜N、P、K施肥總表格產(chǎn)量N施肥量P施肥量K施肥量11.02039137212.072839137214.565639137216.

28、278439137217.7511239137222.5916839137221.6322439137219.3428039137216.1233639137214.113923913726.3922403729.482244937212.462249837214.3322414737217.1022419637221.9422429437222.6422439137221.3422448937222.0722458737224.5322468537215.75224391016.762243914716.892243919316.2422439114017.5622439118619.202

29、2439127917.9722439137215.8422439146520.1122439155819.402243916518.4 GUI人機交互界面程序function varargout = tudougui(varargin)gui_Singleton = 1;gui_State = struct('gui_Name', mfilename, . 'gui_Singleton', gui_Singleton, . 'gui_OpeningFcn', tudougui_OpeningFcn, . 'gui_OutputFcn

30、9;, tudougui_OutputFcn, . 'gui_LayoutFcn', , . 'gui_Callback', );if nargin && ischar(varargin1) gui_State.gui_Callback = str2func(varargin1);end if nargout varargout1:nargout = gui_mainfcn(gui_State, varargin:);else gui_mainfcn(gui_State, varargin:);endfunction tudougui_Openi

31、ngFcn(hObject, eventdata, handles, varargin)handles.output = hObject;guidata(hObject, handles);function varargout = tudougui_OutputFcn(hObject, eventdata, handles) varargout1 = handles.output; function data1_Callback(hObject, eventdata, handles)a=str2num(get(hObject,'String')if(isempty(a) se

32、t(hObject,'String','0')endguidata(hObject,handles);function data1_CreateFcn(hObject, eventdata, handles)if ispc set(hObject,'BackgroundColor','white');else set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor');end function data3_Callba

33、ck(hObject, eventdata, handles)function data3_CreateFcn(hObject, eventdata, handles)q=str2num(get(hObject,'String')if(isempty(a) set(hObject,'String','0')endguidata(hObject,handles);if ispc set(hObject,'BackgroundColor','white');else set(hObject,'Backgroun

34、dColor',get(0,'defaultUicontrolBackgroundColor');endfunction data2_Callback(hObject, eventdata, handles)function data2_CreateFcn(hObject, eventdata, handles)m=str2num(get(hObject,'String')if(isempty(a) set(hObject,'String','0')endguidata(hObject,handles);if ispc set