实验二(插值与拟合)

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,施肥效果分析,- 拟合,数学实验,1,一、实验目的,1.,掌握用MATLAB计算三种插值的方法（编辑Lagr1.m文件,会用拉格朗日程序y=lagr1(x0,y0,x);分段插值程序y=interp1(x0,y0,X)及三次样条插值程序y=interp1(x0,y0,X,spline)或y=spline(x0,y0,X) (其中数组长度自定义(x0,y0同长度,x,y同长度).,2.掌握用MATLAB作线性最小二乘的方法,熟练应用程序语句a=ployfit(x,y,m),x,y是要拟合的数据,是长度自定义的数组,m为拟合多项式的次数,a为输出多项式的系数a=a1,a2,am+!).,3.通过实验学习如何用插值法与拟合法解决实际问题,并注意二者的区别和联系.,2,二、实验内容及要求,1.问题,某地区作物生长所需的主要营养素是氮(N).磷(P).钾(K).现对该地区的土豆与生菜做了一些实验,其数据如下表所示,其中h表示公顷,t表示吨,kg表示公斤.,2.要求,1).试分析施肥量与产量之间的关系,并对所得结果从应用价值与如何改进等方面作出估价.,2).分析实验数据,选择适当的的函数系.在对实测数据进行分析时,应首先根据自己的数据作出散点图(用MATLAB中的函数Plot函数,并分析散点分布情况,选择适当的函数系.,3).运用MATLAB对实测数据进行拟合,并作出拟合函数图;,4).观察并比较误差,提出改进意见或方案.,5).给出MATLAB的程序.,3,表1 土豆产量与施肥量的关系,施肥量,(N)(kg/h),产量(t/h),施肥量,(P)(kg/h),产量(t/h),施肥量,(K)(kg/h),产量(t/h),0,15.18,0,33.46,0,18.98,34,21.36,24,32.47,47,27.35,67,25.72,49,36.06,93,34.86,101,32.29,73,37.96,140,38.52,135,34.03,98,41.04,186,38.44,202,39.45,147,40.09,279,37.73,259,43.15,196,41.26,372,38.43,336,43.46,245,42.17,465,43.87,404,40.83,294,40.36,558,42.77,471,30.75,342,42.73,651,46.22,4,表2 生菜产量与施肥量的关系,施肥量,(N)(kg/h),产量(t/h),施肥量,(P)(kg/h),产量(t/h),施肥量,(K)(kg/h),产量(t/h),0,11.02,0,6.39,0,15.75,28,12.70,49,9.48,47,16.76,56,14.56,98,12.46,93,16.89,84,16.27,147,14.33,140,16.24,112,17.75,196,17.10,186,17.56,168,22.59,294,21.94,279,19.20,224,21.36,391,22.64,372,17.97,280,19.34,489,21.34,465,15.84,336,16.12,587,22.07,558,20.11,392,14.11,686,24.53,651,19.40,5,三、模型假设,1.土豆或生菜生长所需营养素仅为氮,磷,钾;,2.当一个营养素的施肥量变化时,总将另两个营养素的施肥量保持在第7个水平上,如对土豆产量关于磷肥上午施肥量作实验时,氮肥与钾肥的施肥量分别取259kg/h与372kg/h；,3.土豆或生菜的生长不受天气等其他因素的影响；,4.土豆或生菜生长与氮磷钾的关系可以用光滑曲线来近似;,5.当土豆或生菜的产量因为施肥量的增加而大量下降时要停止增加施肥量.,6,四.问题分析,1.引如如下记号,2.根据假设2,当一个营养素的施肥量变化时,总将另两个营养素的施肥量保持在第7个水平上,因此,我们只需考虑某一个营养素的变化与作物产量之间的关系,即利用给出的数据,找出产量与各种营养素之间函数关系.下面,我们只给出土豆产量与施肥量的关系.,7,我们先来研究氮肥与产量的关系,从数据上可以发现:当氮肥增加时产量逐步增加,但当增加到一定时候产量反而减少,这就是农业生产中氮肥的过量使用会造成烧苗的原因.从磷肥与产量的数据可以发现:磷肥的作用很小,当磷肥量为0时,产量为33.46(t/h),而施用(24kg/h)产量反而减少,继续增加施肥量,产量有增有减;最大产量为42.73(t/h),此时的施肥量为342(kg/h),施肥的平均效率为(42.73-33.46)/342=0.0271.钾肥与产量的关系与上述两种肥料有所不同,当施肥量增加到一定的时候,产量增加趋于定值.根据以上分析,我们用曲线回归方法建立土豆产量与施肥量的关系,具体步骤为:,8,1).做出散点图(用MATLAB编程如下):,n=0 34 67 101 135 202 259 336 404 471;,y=15.18 21.36 25.72 32.29 34.03 39.45 43.15 43.46 40.83 30.75;,plot(n,y,r:,n,y,*),xlabel(氮肥的使用量 n),ylabel(土豆的产量 y),title(土豆产量与氮肥使用量的关系),9,10,p=0 24 49 73 98 147 196 245 294 342;,y=33.46 32.47 36.06 37.96 41.04 40.09 41.26 42.17 40.36 42.73;,plot(p,y,c:,n,y,*),xlabel(磷肥的使用量 p),ylabel(土豆的产量 y),title(土豆产量与磷肥使用量的关系),),11,12,k=0 47 93 140 186 279 372 465 558 651;,y=18.98 27.35 34.86 38.52 38.44 37.73 38.43 43.87 42.77 46.22;,plot(p,y,m:,n,y,*),xlabel(钾肥的使用量 k),ylabel(土豆的产量 y),title(土豆产量与钾肥使用量的关系),13,14,从图中可以发现,氮,磷肥的效用可选而次函数来拟合曲线,而钾肥与产量的关系函数可选用指数函数.,利用MATLAB编程(见后附),我们得到产量y与n的函数为:,五. 建立模型,(拟合曲线方程),15,六.模型应用与改进,以上,我们建立了氮.磷.钾肥与产量之间的关系,利用上述函数关系我们可以定量分析计算各种肥料的最佳使用量,我们首先通过市场调查,得到氮磷钾肥和土豆、生菜的市场价格如下表所示：,表3 市场价格（元/吨）,商品,氮肥,磷肥,钾肥,土豆,生菜,价格,350,320,640,800,500,16,由于当一种肥料施肥量改变时，另外的两种肥料都保持在第7个水平上，于是我们可有入下三个方案：,(n,245,465),(259,p,465),(259,245,k),我们只需对上述3个方案求出最大利润,然后比较就可以得出最佳施肥方案:,设每公顷施氮肥n千克的利润为L(n)元,根据题意有,可知,n=327.77时函数有极大值,即最大利润为:L(327.77)=37276.89(元),17,设每公顷施磷肥p千克的利润为L(p)元,根据题意有,最大利润为:L(35.75)=27026.52(元),设每公顷施肥钾k千克的利润为L(k)元,根据题意有,最大利润约为:L(570)=33577.82(元),18,七.结论及建议,综上所述,可知最佳方案为第一个方案(327.77,245,465).,根据农作物的生长原理,氮磷钾3种肥料缺一不可,但又是一个有机整体,因此,要得到农作物的产量与3种肥料之间的使用量的关系,必须考虑3种肥料间的交互影响的数据,也就是说在设计实验时应当采取正交实验,或均匀设计的方法,利用这样实验得到的数据建立的农作物与3种肥料之间的多元函数关系,才能准确地找到最佳施肥方案.,19,附:用MATLAB编写的程序:,n=0 34 67 101 135 202 259 336 404 471;,p=0 24 49 73 98 147 196 245 294 342;,k=0 47 93 140 186 279 372 465 558 651;,yn=15.18 21.36 25.72 32.29 34.03 39.45 43.15 43.46 40.83 30.75;,yp=33.46 32.47 36.06 37.96 41.04 40.09 41.26 42.17 40.36 42.73;,yk=18.98 27.35 34.86 38.52 38.44 37.73 38.43 43.87 42.77 46.22;,plotfit(n,yn,2);,plotfit(p,yp,2);,y=a*(1-b*exp(-c*k);,20,a=41.3:0.1:42.2;b=0.5:0.01:0.59;c=0.09:-0.01:0;,a=41.9:0.1:42.8;b=0.48:0.01:0.57;,c=0.018:-0.001:0.009;,A=a(ones(1,10),:);B=b(ones(1,10),:);C=c(ones(1,10),:);,K=k(ones(1,10),:);,YK=yk(ones(1,10),:);,y=A.*(ones(10,10)-B.*exp(-C.*K);,g=y-YK,w=sum(g.2),Y=inline(153.37+0.64*k-34160*1-0.56*exp(-0.01*k);,%求最小值的程序,x=fminbnd(y,0,651);,ymax=-y(x),21,