数模实验第四版数据拟合与模型参数估计

第1页/共9页数学模型实验—实验报告4学院：河北大学工商学院专业：电气七班姓名：李青青学号：2012484098实验时间：2014/4/15实验地点：B3-301一、实验项目：数据拟合与模型参数估计二、实验目的和要求a.了解数据拟合的原理和Matlab中的有关命令。Polfit：MATLAB函数：p=polyfit(x,y,n)[p,s]=polyfit(x,y,n)说明：x,y为数据点，n为多项式阶数，返回p为幂次从高到低的多项式系数向量p。x必须是单调的。矩阵s用于生成预测值的误差估计。(见下一函数polyval)多项式曲线求值函数：polyval()调用格式：y=polyval(p,x)[y,DELTA]=polyval(p,x,s)说明：y=polyval(p,x)为返回对应自变量x在给定系数P的多项式的值。[y,DELTA]=polyval(p,x,s)使用polyfit函数的选项输出s得出误差估计YDELTA。它假设polyfit函数数据输入的误差是独立正态的，并且方差为常数。则YDELTA将至少包含50%的预测值。Polyval第2页/共9页polyval函数的主要功能是多项式的估值运算，其语法格式为y=polyval(p,x)，输入变量p是长度为n+1的向量，各元素是依次按降幂排列的多项式的系数，函数返回的是那次多项式p在x处的值，x可以是一个数，也可以是一个矩阵或者一个向量,在后两种情况下，该指令计算的是在X中任意元素处的多项式p的估值。polyvalm的主要功能是用于matlab中多项式求值。其语法格式为y=polyvalm（a，A），其中a为多项式行向量表示，A为指定矩阵。Lsqlin约束线性最小二乘函数lsqlin格式x=lsqlin(C,d,A,b)%求在约束条件下，方程Cx=d的最小二乘解x。x=lsqlin(C,d,A,b,Aeq,beq)%Aeq、beq满足等式约束，若没有不等式约束，则设A=[]，b=[]。x=lsqlin(C,d,A,b,Aeq,beq,lb,ub)%lb、ub满足，若没有等式约束，则Aeq=[]，beq=[]。x=lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0)%x0为初始解向量，若x没有界，则lb=[]，ub=[]。x=lsqlin(C,d,A,b,Aeq,beq,lb,ub,x0,options)%options为指定优化参数lsqcurvefit第3页/共9页最常见的调用格式如下：X=LSQCURVEFIT(FUN,X0,XDATA,YDATA).其中FUN为一个函数，已M文件或匿名函数存在。若FUN以M文件形式存在，那么FUN在调用语句中的格式为：@(x,xdata)FUN(x,xdata,c).@(x,xdata)中分别表示待求参数，xdata表示的是自变量,c是可以传递到函数里面的常数。Lsqnonlin：lsqnonlin解决非线性最小二乘问题，包括非线性数据拟合问题而不是计算的值f（x）（平方和），需要用户定义函数lsqnonlin求向量值函数然后，在矢量的术语，你可以重申这一优化问题其中x是一个向量和f（x）是一个函数，返回一个向量值。X=lsqnonlin（乐趣，X0）开始在点X0并找到一个最小的有趣的功能描述的平方和。快乐应该返回一个向量值不值的平方和。（算法隐含和广场的乐趣（X）。）X=lsqnonlin（乐趣，X0，LB，UB）定义了一组上下对X设计变量的范围，所以，解决方案总是在范围≤x≤UB的LB。X=lsqnonlin（乐趣，X0，LB，UB，选项）最大限度地减少结构中的优化选项指定的选项。使用optimset设置这些选项。通过空矩阵的LB和UB如果没有界限的存在。【X，resnorm]=lsqnonlin（……）返回x的平方范数的剩余价值的总和（意思：（X）。^2）。第4页/共9页【X，resnorm，残余]=lsqnonlin（……）返回剩余的娱乐价值观（x）在解xx=lsqnonnegx=lsqnonneg（C，D）返回向量x的最小范数（C*XD）受x=0。C和D必须是真实的。x=lsqnonneg（C，D，x0X0）使用为出发点，如果所有的X0=0；否则，则使用默认值。默认的出发点是原点（默认是用来当X0==[]或只有两个输入参数提供）。x=lsqnonneg（C，D，X0，选项）最大限度地减少在结构优化参数指定的选项。你可以使用optimset函数定义这些参数。lsqnonneg使用这些选项的结构域：displayLevel显示。”“不显示输出；最后的“仅显示最终的输出；“通知”（默认）显示输出只有当函数不converge.tolxtermination公差对X【X，resnorm]=lsqnonneg（……）返回的平方范数的剩余价值：规范（C*XD）^2b.练习模型参数估计三、实验内容根据美国人口从1790年到1990年间的人口数据（如下表），确定人口指数增长模型和Logistic模型中的待定参数，估计出美国2010年的人口，同时第5页/共9页画出拟合效果的图形。表1美国人口统计数据年份1790180018101820183018401850人口(×106)3.95.37.29.612.917.123.2年份1860187018801890190019101920人口(×106)31.438.650.262.976.092.0106.5年份193019401950196019701980人口(×106)123.2131.7150.7179.3204.0226.5提示：指数增长模型：rtextx0)(Logistic模型：011mrtmxxtxex解：模型一:指数增长模型。Malthus模型的基本假设下，人口的增长率为常数，记为r，记时刻t的人口为)(tx，（即)(tx为模型的状态变量）且初始时刻的人口为0x，因为0)0(xxrxdtdx由假设可知0()rtxtxe经拟合得到：第6页/共9页2120010120()ln()ln,ln(),,lnrtayataxtxextxrtraxeyxtarax程序：t=1790:10:1980;x(t)=[3.95.37.29.612.917.123.231.438.650.262.976.092.0106.5123.2131.7150.7179.3204.0226.5];y=log(x(t));a=polyfit(t,y,1)r=a(1),x0=exp(a(2))x1=x0.*exp(r.*t);plot(t,x(t),'r',t,x1,'b')结果：a=0.0214-36.6198r=0.0214x0=1.2480e-016所以得到人口关于时间的函数为：0.02140()txtxe，其中x0=1.2480e-016，输入：t=2010;x0=1.2480e-016;x(t)=x0*exp(0.0214*t)得到x(t)=598.3529。即在此模型下到2010年人口大约为598.3529610。第7页/共9页17801800182018401860188019001920194019601980050100150200250300350模型二：阻滞增长模型（或Logistic模型）由于资源、环境等因素对人口增长的阻滞作用，人口增长到一定数量后，增长率会下降，假设人口的增长率为x的减函数，如设)/1()(mxxrxr，其中r为固有增长率(x很小时)，mx为人口容量（资源、环境能容纳的最大数量），于是得到如下微分方程：0)0()1(xxxxrxdtdxm建立函数文件curvefit_fun2.mfunctionf=curvefit_fun2(a,t)f=a(1)./(1+(a(1)/3.9-1)*exp(-a(2)*(t-1790)));在命令文件main.m中调用函数文件curvefit_fun2.m%定义向量（数组）x=1790:10:1990;y=[3.95.37.29.612.917.123.231.438.650.262.976...第8页/共9页92106.5123.2131.7150.7179.3204226.5251.4];plot(x,y,'*',x,y);%画点，并且画一直线把各点连起来holdon;a0=[0.001,1];%初值%最重要的函数，第1个参数是函数名（一个同名的m文件定义），第2个参数是初值，第3、4个参数是已知数据点a=lsqcurvefit('curvefit_fun2',a0,x,y);disp(['a='num2str(a)]);%显示结果%画图检验结果xi=1790:5:2020;yi=curvefit_fun2(a,xi);plot(xi,yi,'r');%预测2010年的数据x1=2010;y1=curvefit_fun2(a,x1)holdoff运行结果：a=311.95310.02798178y1=267.1947其中a(1)、a(2)分别表示011mrtmxxtxex中的mx和r，y1则是对美国美国2010年的人口的估计。第9页/共9页1750180018501900195020002050050100150200250300