插值拟合方法

第4讲数据的插值拟合法4.1数据的拟合方法4.2数据的插值方法4.1数据的拟合方法已知n个数据点（xi,yi)，i=1,2,…,n，xi互不相同，如何寻求函数y=f(x)，使f(x)在某种准则下与这n个点最接近？

拟合模型通过寻找简单的因果变量之间的数量关系，对未知的情形作出预测与预报.Remark：不必要求近似函数的曲线或曲面通过所有的数据点.一、拟合原理问题二、拟合模型的分类与方法1.直线拟合：用一次函数或称线性函数拟合数据.2.曲线拟合：若直线拟合效果不佳，可提高拟合精度，用曲线拟合数据.常用的是二次函数、三次函数等高次多项式，有时也会用到指数函数、对数函数、三角函数等.3.观察数据修匀：根据数据分布的总趋势去剔除观察数据中的偶然误差，即数据修匀（或称数据光滑）问题.4.分段拟合：在不同段上用不同的低次多项式进行拟合.选择何种曲线拟合最好呢？首先可在坐标轴上画出数据的散点图，通过观察选择几种合适的曲线分别拟合，通过比较，哪条曲线的最小二乘指标J最小即为最好的拟合曲线.表示第i个点的拟合值与实际值的绝对误差第i个点的拟合值第i个点的实际值一、模型假设假设所测数据均为抽样数据.二、模型的分析与建立先画出散点图，见下图.问题一

【温度与电阻的关系模型】有一个对温度敏感的电阻，现测得一组温度t与电阻R的数据。见表9.1t20.532.751.07395.7R7658268739421032试给出温度与电阻间的函数关系，并计算温度为60度时的电阻值表9.1t=[20.532.7517395.7];R=[7658268739421032];plot(t,R,'*')xlabel('t')ylabel('R')观察上图不难发现：散点基本上在一条直线上，因此可假设电阻与温度满足一次函数（或称线性函数）.

设拟合函数为y=a1+a2xt=[20.532.7517395.7];R=[7658268739421032];A=[ones(5,1),t'];b=R';a1a2=A\ba1a2=702.09683.3987三、模型求解可用Matlab计算.解法一：使用regress()函数t=[20.532.7517395.7];R=[7658268739421032];x=[ones(5,1),t'];y=R';b=regress(y,x,0.05)b=702.09683.3987解法二：使用ployfit()函数

cleart=[20.532.7517395.7];R=[7658268739421032];x=t';y=R';p=polyfit(x,y,1)p=3.3987702.0968即a1=702.0968,a2=3.3987，因此拟合的函数为y=702.0968+3.3987x将x=60代入计算得y=906cleart=[20.532.7517395.7];R=[7658268739421032];plot(t,R,'*')xlabel('t');ylabel('R')holdont=20:0.1:100;plot(t,702.0968+3.3987.*t)实际值与拟合图regress()

函数主要用于线性拟合，在拟合时进行显著性检验，故称为回归函数.Polyfit()函数主要是利用多项式拟合.它可以是线性或非线性.Remark:polyfit(x,y,m)表示用m次多项式拟合数据x,y.

regress()

函数与polyfit()

函数的区别问题2

【农业生产实验模型】在研究农业生产的试验中，为分析某地区土豆产量与化肥的关系，得到了每公顷地的氮肥的施肥量与土豆产量的对应关系，见图9.2氮肥量(kg)03467101135202259336404471土豆产量(kg)15.1821.3625.7232.2934.0339.4543.1543.4640.8330.75请根据表9.2的数据，给出土豆产量与氮肥施肥量之间的关系表9.2一、模型假设

假设实验数据为抽样数据.

假设其他化肥用量不变.二、模型分析与建立先用Matlab画出散点图.clearx=[03467101135202259336404471];y=[15.1821.3625.7232.2934.0339.4543.1543.4640.8330.75];plot(x,y,'+')xlabel('x(danfeiliang)')ylabel('y(tudouchanliang)')可以看出散点图呈二次曲线图形，可取如下拟合函数其中x表示氮肥量，y表示土豆产量。a,b,c为待定系数。三、模型求解abc=-0.00030.197114.7416clearx=[03467101135202259336404471];y=[15.1821.3625.7232.2934.0339.4543.1543.4640.8330.75];abc=polyfit(x,y,2)x1=0:471;y1=polyval(abc,x1);plot(x,y,'*',x1,y1)xlabel('x(danfeiliang)')ylabel('y(tudouchanliang)')问题3

【血药浓度模型】通过实验测得一次性快速静脉注射300mg药物后的血药浓度数据，见表9.3t(h)0.250.511.523468y(ug/ml)19.2118.1515.3614.1012.899.327.455.243.01求血药深度随时间的变化规律y(t)表9.3一、模型假设1.假设t=0时，y=0.2.假设实验数据为抽样数据，能反映血药浓度与时间的关系二、模型分析与建立

利用Matlab画出散点图.cleart=[0.250.511.523468];y=[19.2118.1515.3614.1012.89 9.327.455.243.01];plot(t,y,'o')xlabel('t(shijian)')ylabel('y(nongdu)')由图可见，散点图大致呈负指数函数形态，可令其中a,b>0为待定系数.三、模型求解解法一：将其线性化处理.将等式

两边同时取对数，得令则方程变为用Matlab求解如下:cleart=[0.250.511.523468];y=[19.2118.1515.3614.1012.89 9.327.455.243.01];Y=log(y);b=polyfit(t,Y,1)b=-0.23472.9943即b1=2.9943b2=-0.2347则血药浓度与时间的关系为cleart=[0.250.511.523468];y=[19.2118.1515.3614.1012.89 9.327.455.243.01];Y=log(y);b=polyfit(t,Y,1);a=exp(b(2));b=-b(1);plot(t,y,'*');holdont=0:0.01:9;plot(t,a.*exp(-b.*t),'r')xlabel('t')ylabel('y')解法2

使用非线性拟合函数nlinfit().首先建立Fun7_1.m文件如下：functiony=Fun7_1(a,t)y=a(1)*exp(-a(2)*t);不妨取a=20,由t=1时，y=15.36算出b=0.264取a0=[200.264];用Matlab求解如下cleart=[0.250.511.523468];y=[19.2118.1515.3614.1012.89 9.327.455.243.01];a0=[20,0.264];a=nlinfit(t,y,'fun7_1',a0)a=20.24130.2420t=[0.250.511.523468];y=[19.2118.1515.3614.1012.89 9.327.455.243.01];plot(t,y,'*');holdona0=[20,0.264];a=nlinfit(t,y,'fun7_1',a0);t=0:0.01:9;plot(t,a(1).*exp(-a(2).*t))holdonplot(t,19.9714.*exp(-0.2347.*t),'r')xlabel('t');ylabel('y')nlinfit函数的用法

nlinfit函数采用的迭代法，其中a0为迭代初值.若所给初值离最优解比较近，则迭代求出该最优解的概率就很高.

如何估计初值，暂无确切方法.可以在得到解后画出函数图形，看看实验点是否都在曲线附近.若相差太大，可以考虑重新给出初值重新计算.问题4一、模型假设假设氯气与生产时间之间满足其中a,b为待定系数.二、模型分析与建立此问题实质上是确定待定系数a,b的值.三、模型求解首先定义非线性函数fun7_2.m文件：functiony=fun7_2(beta0,x)a=beta0(1);b=beta0(2);y=a+(0.49-a)*exp(-b*(x-8));然后在命令窗口中输入clearx=[88101010101212121214141416161618182020202022222424242626262828303030323234363638384042];y=[0.490.490.480.470.480.470.460.460.450.430.450.430.430.440.430.430.460.450.420.420.430.410.410.40.420.40.40.410.40.410.410.40.40.40.380.410.40.40.410.380.40.40.390.39];beta0=[0.300.20];ab=nlinfit(x,y,'fun7_2',beta0)ab=0.39040.1028即a=0.3904,b=0.1028.所以模型函数为

作业：用Matlab画出相关图形.问题5一、模型假设假设中国人口的变化满足一定规律.二、模型分析、建立与求解首先画出散点图.clearx=1949:5:1994;y=[5.466.778.19.19.810.311.311.8];plot(x,y,'*')xlabel('x(nianfen)')ylabel('y(renkoushu)')三、模型求解模型一：用线性函数拟合.设我国人口数量满足以下模型y=a+bx其中a,b为待定系数.通过计算得y=-283.232+0.148x模型二：用指数函数拟合.设我国人口数量满足以下模型y=aebx其中a,b为待定参数.模型两边取对数得lny=lna+bxy=4.1444*10-15e0.0179xx=1949:5:1994;y=[5.466.778.19.19.810.311.311.8];a=polyfit(x,y,1);x1=1949:0.1:1994;y1=a(2)+a(1)*x1;b=polyfit(x,log(y),1);y2=exp(b(2))*exp(b(1)*x1);plot(x,y,'*')holdonplot(x1,y1)holdonplot(x1,y2,'r')xlabel('x(nianfen)')ylabel('y(renkoushu)')用两个模型分别计算相应年度的人口数及其误差见下表.分别计算两个模型的最小二乘指标得从而线性模型更适合中国人口的增长.4.2数据的插值方法若知道函数y=f(x)在n个互异的点x1,x2,…,xn的函数值y1,y2,…,yn

.如何估计此函数在另一点a的函数值？一、插值原理问题考虑构造一个过x1,x2,…,xn的次数不超过n的多项式y=Ln(x)，使其满足Ln(xk)=yk,k=1,2,…,n

然后用Ln(a)作为准确值L(a)的近似值.这种方法叫作插值.

Remark:插值方法要求近似函数经过所有的已知点.插值方法一般有：拉格朗日（Lagrange）插值、分段线性插值、埃尔米特（Hermite）插值、样条插值以及分形插值等.Matlab的插值函数分为内部插值与外部插值.内部插值：要求已知点x是单调的，且被插值点xi不能超过x的范围.如interp1()、interp2()、interpn()等.而griddata()既可以计算内部插值，也可以计算外部插值.二、插值方法数据拟合要求得到一个具体的近似函数表达式，而数据插值不一定得到近似函数的表达式，仅通过插值方法找到未知点对应的近似函数值.数据拟合得到的近似函数表达式不一定要经过所有已知点，而数据插值要求近似函数经过已知的所有数据点.数据拟合与数据插值的区别一、模型假设假设飞机机翼截面下轮廓的变化是连续的.二、模型分析与建立画出散点图问题7从散点图很难观察出函数形式,用插值的方法.由于只有一个变量且是单调的，可用一元插值函数interp1()进行计算.三、模型求解用Matlab求解如下：x1=[035791112131415];y1=[01.21.722.121.81.211.6];x=0:0.1:15;y=interp1(x1,y1,x,'spline')plot(x1,y1,'.',x,y,'r')gridtitle('spline')xlabel('x')ylabel('y')插值后的结果如下图显示.一、模型假设假设电压的变化是连续的.二、模型分析与建立画出散点图.问题8由散点图可观察出电容器两端的电压随充电时间的关系大致呈对数关系，本题不要求给出函数形式，因此同样使用插值的方法.三、模型求解用Matlab求解如下：插值后结果如图所示.functionwenti8_2x1=[12346.5912];y1=[6.27.38.299.610.110.4];x=0:0.1:12;y=interp1(x1,y1,x,'spline')plot(x1,y1,'.',x,y,'r')gridtitle('spline')xlabel('x')ylabel('y')Remark:一元插值函数interp1的基本调用格式为interp1(x,y,cx,‘method’).

其中x,y分别表示已知数据点的横、纵坐标.x必须单调，cx为需要插值的横坐标，method为可选参数，可为以下四个值之一：

1）nearest——最近邻点插值

2）linear——线性插值（可缺省）

3）spline——三次样条插值

4）cubic——三次插值一、模型假设假设温度的变化是连续的.二、模型分析与建立首先画出散点图.问题9从散点图难以观察出函数形式,下面用插值的方法.三、模型求解用Matlab求解如下.x1=1:1:12;y1=[589152529313022252724];x=[3.26.57.111.7];y=interp1(x1,y1,x)plot(x1,y1,'*',x,y,'O')xlabel('x')ylabel('y')y=10.200030.000030.900024.9000运行结果如图所示.问题10一、模型假设1.假设河床的深度是连续变化的.2.假设紧靠河床铺设电缆，即电缆长度等于河床长度.二、模型分析与建立画出河床观测的散点图x=5:5:100;y=[2.412.962.152.653.124.235.126.215.68...4.223.913.262.852.353.023.634.123.462.080];y1=10-y;plot(x,y1,'*')axis([0100010])xlabel('x')ylabel('y')gridon

三、模型求解利用分段线性插值，并可以在此基础上利用梯形法求积分命令trapz计算河床面积.同时利用每段连续线长度之和来近似河床曲线长度.x=5:5:100;y=[2.412.962.152.653.124.235.126.215.68...4.223.913.262.852.353.023.634.123.