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

数学模型实验—实验报告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得出误差估计Y DELTA。

它假设polyfit函数数据输入的误差是独立正态的，并且方差为常数。

则Y DELTA将至少包含50%的预测值。

Polyvalpolyval函数的主要功能是多项式的估值运算，其语法格式为y = poly val(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最常见的调用格式如下：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）。

【X，resnorm，残余] = lsqnonlin（……）返回剩余的娱乐价值观（x）在解xx = lsqnonnegx = lsqnonneg（C，D）返回向量x的最小范数（C * XD）受x > = 0。

C 和D必须是真实的。

x = lsqnonneg（C，D，x0 X0）使用为出发点，如果所有的X0 > = 0；否则，则使用默认值。

默认的出发点是原点（默认是用来当X0 = = []或只有两个输入参数提供）。

x = lsqnonneg（C，D，X0，选项）最大限度地减少在结构优化参数指定的选项。

你可以使用optimset函数定义这些参数。

lsqnonneg使用这些选项的结构域：displayLevel显示。

”“不显示输出；最后的“仅显示最终的输出；“通知”（默认）显示输出只有当函数不converge.tolxtermination公差对X【X，resnorm ] = lsqnonneg（……）返回的平方范数的剩余价值：规范（C * XD）^ 2b. 练习模型参数估计三、实验内容根据美国人口从1790年到1990年间的人口数据（如下表），确定人口指数增长模型和Logistic模型中的待定参数，估计出美国2010年的人口，同时画出拟合效果的图形。

表1 美国人口统计数据提示：指数增长模型：rte x t x 0)(=Logistic 模型：()011mrtm x x t x e x -=⎛⎫+- ⎪⎝⎭解：模型一:指数增长模型。

Malthus 模型的基本假设下，人口的增长率为常数，记为r ，记时刻t 的人口为 )(t x ，（即)(t x 为模型的状态变量）且初始时刻的人口为0x ，因为⎪⎩⎪⎨⎧==0)0(x x rxdt dx 由假设可知0()rt x t x e = 经拟合得到：}2120010120()ln ()ln ,ln (),,ln rt a y a t a x t x e x t x rt r a x ey x t a r a x =+=⇒=+⇒===== 程序：t=1790:10:1980;x(t)=[3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76.0 92.0 106.5123.2 131.7 150.7 179.3 204.0 226.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.0214 x0= 1.2480e-016所以得到人口关于时间的函数为：0.02140()t x t x e =，其中x0 = 1.2480e-016， 输入：t=2010;x0 = 1.2480e-016; x(t)=x0*exp(0.0214*t)得到x(t)= 598.3529。

即在此模型下到2010年人口大约为598.3529 610⨯。

17801800182018401860188019001920194019601980050100150200250300350模型二：阻滞增长模型（或 Logistic 模型） 由于资源、环境等因素对人口增长的阻滞作用，人口增长到一定数量后，增长率会下降，假设人口的增长率为 x 的减函数，如设)/1()(m x x r x r -=，其中 r 为固有增长率 (x 很小时 ) ，m x 为人口容量（资源、环境能容纳的最大数量）， 于是得到如下微分方程： ⎪⎩⎪⎨⎧=-=0)0()1(xx x x rx dtdxm 建立函数文件curvefit_fun2.m function f=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.9 5.3 7.2 9.6 12.9 17.1 23.2 31.4 38.6 50.2 62.9 76 ...92 106.5 123.2 131.7 150.7 179.3 204 226.5 251.4]; plot(x,y,'*',x,y); % 画点，并且画一直线把各点连起来 hold on;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) hold off 运行结果：a=311.9531 0.02798178 y1 =267.1947其中a(1)、a(2)分别表示()011mrtm x x t x e x -=⎛⎫+- ⎪⎝⎭中的m x 和r ，y1则是对美国美国2010年的人口的估计。

1750180018501900195020002050050100150200250300。