数学实验8-曲线拟合与插值

MathematicsLaboratory阮小娥博士ExperimentsinMathematics数学实验在实际问题中，我们常常会遇到下列问题（1）变量间存在函数关系，但只能给出一离散点列上的值.例如，从实验中得到一个数据表,或是一组观测数据.（2）变量间的函数关系可以表示,但计算复杂,只能计算特殊点的函数值.例如，求指数函数、对数函数、三角函数、反三角函数值等.为了研究自变量与因变量间的变化关系，我们需要建立变量间的函数关系，从而可以计算原始数据以外需要处的值.解决这类问题的方法：数据拟合、数据插值实验13人口数量预测模型实验2、掌握在最小二乘意义下数据拟合的理论和方法.1、学会用MATLAB软件进行数据拟合3、通过对实际问题的分析和研究，初步掌握建立数据拟合数学模型的方法实验目的据人口统计年鉴，知我国从1949年至1994年人口数据资料如下：(人口数单位为：百万)（1）在直角坐标系上作出人口数的图象。（2）建立人口数与年份的函数关系，并估算1999年的人口数。实验问题年份

1949

1954

1959

1964

1969

人口数

541.67

602.66

672.09704.99

806.71

年份1974

1979

1984

1989

1994人口数

908.59

975.421034.75

1106.76

1176.74

如何确定a,b?线性模型1曲线拟合问题的提法:

已知一组（二维）数据，即平面上的n个点),(iiyx，

ixni,,,2,1L=互不相同，寻求一个函数（曲线）)(xfy=，使)(xf在某种准则下与所有数据点最为接近，即曲线拟合得最好，如图:

xy0++++++++一、曲线拟合确定f(x)使得

达到最小

最小二乘准则

２.用什么样的曲线拟合已知数据?常用的曲线函数系类型：(１)画图观察(２)理论分析指数曲线：

双曲线（一支):

多项式：

直线：

３拟合函数组中系数的确定二、人口预测的线性模型

对于开始提出的实验问题,代如数据，计算得从而得到人口数与年份的函数关系为把x=1999代如，估算出1999年的人口数为y=1252.1（百万）＝12.52亿1999年实际人口数量为１２.６亿。线性预测模型英国人口学家Malthus根据百余年的人口统计资料，于1798年提出了著名的人口自然增长的指数增长模型。三、人口预测的Malthus模型基本假设

:人口(相对)增长率r

是常数x(t)~时刻t的人口,t=0时人口数为x0指数增长模型实际中，常用1.由前100年的数据求出美国的人口增长Malthus模型。2.预测后100年（每隔10年）的人口状况。3.根据预测的人口状况和实际的人口数量,讨论人口模型的改进情况。美国1790年－1980年每隔10年的人口记录226.5204.0179.3150.7131.7123.2106.592.076.062.9人口(百万)1980197019601950194019301920191019001890年份50.238.631.423.217.112.99.67.25.33.9人口(百万)1880187018601850184018301820181018001790年份例１解：取得最小值.其中,表示人口数量。表示年份,解方程组:即得参数的值.使得问题转化为求参数％prog41.m%%Thisprogramistopredictthenumberofpopulation%formatlongt1=[1790;1800;1810;1820;1830;1840;1850;1860;1870;1880];t2=[1890;1900;1910;1920;1930;1940;1950;1960;1970;1980];x1=[3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2];x2=[62.9;76.0;92.0;106.5;123.2;131.7;150.7;179.3;204.0;226.5];lnx1=log(x1);lnx2=log(x2);a12=sum(t1);a11=10;a21=a12;a22=sum(t1.^2);d1=sum(lnx1);d2=sum(lnx1.*t1);

A=[a11,a12;a21,a22];D=[d1;d2];

ab=inv(A)*D;

disp('a=');disp(ab(1));

disp('b=');disp(ab(2));

fori=1:10

xx1(i)=exp(ab(1)+ab(2)*t1(i));

end

fori=1:10

xx2(i)=exp(ab(1)+ab(2)*t2(i));

end

plot(t1,x1,'r*--',t1,xx1,'b+-',t2,x2,'g*--',t2,xx2,'m+-');

a=-49.79535457790735b=0.02859807120038仿真结果表明：人口增加的指数模型在短期内基本上能比较准确地反映人口自然增长的规律，但长期预测误差很大，需要修正预测模型。拟合曲线原始数据曲线四、人口预测的Logistic模型人口增长到一定数量后，增长率下降的原因：资源、环境等因素对人口增长的阻滞作用且阻滞作用随人口数量增加而变大假设r~固有增长率(x很小时)k~人口容量（资源、环境能容纳的最大数量）r是x的减函数例１的Logistic模型留给同学们练习五、多项式拟合的Matlab指令a=polyfit（xdata，ydata，n）其中n表示多项式的最高阶数

xdata，ydata为要拟合的数据，它是用向量的方式输入。输出参数a为拟合多项式

y=a1xn+…+anx+an+1的系数a=[a1,…,an,an+1]。多项式在x处的值y可用下面程序计算。

y=polyval(a,x)

用多项式拟合人口模型%Thisprogramistopredictthemodelofpopulationby4-degreepolynomial%%prog42.m%formatlongt1=[1790;1800;1810;1820;1830;1840;1850;1860;1870;1880];t2=[1890;1900;1910;1920;1930;1940;1950;1960;1970;1980];t=[t1;t2];P1=[3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2];P2=[62.9;76.0;92.0;106.5;123.2;131.7;150.7;179.3;204.0;226.5];P=[P1;P2];n=4;%Thedegreeofthefittingpolynomial%[a,s]=polyfit(t1,P1,n);y=polyval(a,t);%aisthecoefficientsvectorfromn-degreeto0-degree%plot(t,P,'r*--',t,y,'b+-');23a=1.0e+006*-0.000000000000140.00000000107892-0.000003048785950.00381927346813-1.79012132225427仿真结果表明,人口增加的模型用多项式拟合能比较准确地反映人口自然增长的规律，对长期预测具有指导意义。例2:海底光缆线长度预测模型某一通信公司在一次施工中,需要在水面宽为20m的河沟底沿直线走向铺设一条沟底光缆.在铺设光缆之前需要对沟底的地形做初B2468101214161820986420ADC探测到一组等分点位置的深度数据如下表所示.25步探测,从而估计所需光缆的长度,为工程预算提供依据.基本情况如图所示.10.9310.809.818.867.957.959.1510.2211.2912.6113.32201918171615141312111013.2812.2611.1810.139.058.027.967.968.969.01深度(m)9876543210分点21个等分点处的深度(1)预测通过这条河沟所需光缆长度的近似值.(2)作出铺设沟底光缆的曲线图.解：用12次多项式函数拟合光缆走势的曲线图如下仿真结果表明,拟合曲线能较准确地反映光缆的走势图.Thelengthofthelabelis

L=26.3809(m)假设所铺设的光缆足够柔软,在铺设过程中光缆触地走势光滑,紧贴地面,并且忽略水流对光缆的冲击.%prog45.mThisprogramistofitthedatabypolynomial%formatlongt=linspace(0,20,21);x=linspace(0,20,100);P=[9.01,8.96,7.96,7.97,8.02,9.05,10.13,11.18,12.26,13.28,13.32,12.61,11.29,10.22,9.15,7.90,7.95,8.86,9.81,10.80,10.93];[a,s]=polyfit(t,P,12);yy=polyval(a,x);plot(x,yy,'r*--',t,P,'b+-');L=0;fori=2:100L=L+sqrt((x(i)-x(i-1))^2+(yy(i)-yy(i-1))^2);enddisp('ThelengthofthelabelisL=');disp(L);formatlongt=linspace(0,20,21);x=linspace(0,20,100);P=[9.01,8.96,7.96,7.97,8.02,9.05,10.13,11.18,12.26,13.28,13.32,12.61,11.29,10.22,9.15,7.90,7.95,8.86,9.81,10.80,10.93];n=input(‘n=’)%通过键盘输入拟合次数[a,s]=polyfit(t,P,n);yy=polyval(a,x);p1=polyval(a,t);d=norm(P-p1)%计算拟合误差plot(x,yy,'r*--',t,P,'b+-');L=0;fori=2:100L=L+sqrt((x(i)-x(i-1))^2+(yy(i)-yy(i-1))^2);enddisp('ThelengthofthelabelisL=');disp(L);六、最小二乘曲线拟合有一组数据xi，yi（i=1,2,…,n）满足某一函数原型，其中a为待定系数向量，求出一组待定系数的值使得目标函数最小：最小二乘曲线拟合函数lsqcurvefit的调用格式：[a,Jm]=lsqcurvefit(Fun,a0,x,y)Fun为原型函数的matlab表示，可以是M-函数或inline()函数a0为最优化初值x和y为原始输入输出数据向量a为返回的待定系数向量Jm为在待定系数下目标函数的值例3

已知数据可能满足求满足数据的最小二乘解a、b、c和d的值.x0.10.20.30.40.5y2.32012.64702.97073.28853.6008x0.60.70.80.91.0y3.90904.21474.51914.82325.1275输入已知数据点：x=0.1:

0.1

:

1;y=[2.3201,2.6470,2.9707,3.2885,3.6008,3.9090,4.2147,4.5191,4.8232,5.1275];编写函数functiony=f3(a,x)y=a(1)*x+a(2)*x.^2.*exp(-a(3)*x)+a(4);待定系数求解a=lsqcurvefit('f3',[1;2;2;3],x,y);a'ans=2.45752.45571.44372.0720绘制曲线：>>y1=f3(a,x);plot(x,y,x,y1,’o’)插值问题：实验14插值问题插值条件插值函数插值节点如果是多项式,则称为插值多项式求插值函数的方法称为插法…………[a,b]称为插值区间如何构造P(x)?1、多项式插值法当n=0时,只有一个插值节点的情形当n=1时,有两个插值节点的情形当n=2时,有三个插值节点的情形是否任意给定n+1个不同的插值节点都可以构造出满足插值条件的插值多项式?由克莱姆法则知方程组有唯一解，即满足(1.1)的插值多项式存在且唯一。插值多项式的存在唯一定理：在次数不超过n的多项式集合中，满足插值条件的插值多项式是存在并且唯一的。拉格朗日(Lagrange)插值多项式表示插值多项式的最紧凑的方式是拉格朗日形式：

Lagrange插值多项式的优点在于不要求数据点是等间隔的，其缺点是数据点数不宜过大，通常不超过7个，否则计算工作量大且误差大，计算不稳定。分段线性插值分段线性插值示意图…………分段二次插值示意图…………分段二次插值三次样条插值函数定义对于区间[a,b]上给定的一个分划如果函数S(x)在子区间上都是不超过3次的多项式，并且2阶导数在内节点处连续，则称为区间[a,b]上以为节点的三次样条函数。对于函数，若还满足插值条件：则称为在区间上的三次样条插值函数。是在飞机或轮船等的设计制造过程中为描绘出光滑的外形曲线(放样)所用的工具.2、插值函数的MATLAB实现（1）interp1函数interp1函数的调用格式：y=interp1(x0,y0,x,method)其中:1）x0、y0为样本点，y为插值点自变量值x对应的函数值。（2）method共有6种参数可供选择，当省略method时，即默认为linear线性插值。线性插值：method=‘linear’三次样条插值:method=‘spline’;立方插值：method=‘pchip’or‘cubic’例4已知某转子流量计在100～1000mL/min流量范围内,刻度值与校正值的关系如表所示.试用线性插值法计算流量计的刻度值为785时,实际流量为多少？刻度值校正值刻度值校正值100105.3600605.8200207.2700707.4300308.1800806.7400406.9900908.0500507.51000107.9解：Matlab计算程序如下：X=[100,200,300,400,500,600,700,800,900,1000];Y=[105.3,207.2,308.1,406.9,507.5,605.8,707.4,806.7,908.0,107.9];Xk=780;Yk=interp1(X,Y,Xk)执行结果：Yk=786.8400这里：X和Y分别表示样本点的刻度值和校正值；

Xk和Yk分别表示插值点的刻度值和校正值。功能

三次样条数据插值格式

y=spline(x0,y0,x)

与y=interp1(x0,y0,x,’spline’)等价，其中参数x0、y0、x，y的意义及要求与线性插值interp1中的完全一致。（2）spline函数例5对离散分布在y=exp(x)sin(x)函数曲线上的数据点进行样条插值计算：

x=[024581212.817.219.920];y=exp(x).*sin(x);xx=0:0.25:20;yy=spline(x,y,xx);plot(x,y,'bo',xx,yy,'r*')例6

已知某型号飞机的机翼断面下缘轮廓线上的部分数据如表所示:假设需要得到x

坐标每改变0.1时的y

坐标,

分别用两种插值方法对机翼断面下缘轮廓线上的部分数据加细,并作出插值函数的图形.xy0031.251.772.092.1112.0121.8131.2141.0151.6程序如下：x=[0,3,5,7,9,11,12,13,14,15]y=[0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6]xi=[0:0.1:15]yi=interp1(x,y,xi,'spline')zi=spline(x,y,xi)plot(xi,yi,’bO’,xi,zi,’r*’)（3）lagrange插值函数function

y=lagrange(x0,y0,x)ii=1:length(x0);y=zeros(size(x));fori=iiij=find(ii~=i);y1=1;forj=1:length(ij),y1=y1.*(x-x0(ij(j)));endy=y+y1*y0(i)/prod(x0(i)-x0(ij));end例7对进行Lagrange插值x0=-1+2*[0:10]/10;y0=1./(1+25*x0.^2);x=-1:0.01:1;y=lagrange(x0,y0,x);ya=1./(1+25*x.^2);

plot(x,ya,'*',x,y,'O')3、二维网格数据的插值问题（1）二维插值函数interp2的调用格式：zi=interp2(x0,y0,z0,xi,yi)zi=interp2(x0,y0,z0,xi,yi,

method)*临近点插值：method=‘nearest’*线性插值：method=‘linear’（缺省算法）*三次样条插值：method=‘spline’*立方插值：method=‘pchip’or‘cubic’例8

由二元函数获得一些较稀疏的网格数据，对整个函数曲面进行各种插值，并比较插值结果%绘制已知数据的网格图[x,y]=meshgrid(-3:0.6:3,-2:0.4:2);z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);surf(x,y,z);axis([-3,3,-2,2,-0.7,1.5])%默认线性插值算法进行插值[x1,y1]=meshgrid(-3:.2:3,-2:.2:2);z1=interp2(x,y,z,x1,y1);surf(x1,y1,z1),axis([-3,3,-2,2,-0.7,1.5])%立方和样条插值z2=interp2(x,y,z,x1,y1,'cubic');z3=interp2(x,y,z,x1,y1,'spline');surf(x1,y1,z2),axis([-3,3,-2,2,-0.7,1.5])figure;surf(x1,y1,z3),axis([-3,3,-2,2,-0.7,1.5])[x,y]=meshgrid(-3:0.6:3,-2:0.4:2);z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);subplot(2,2,1);surf(x,y,z);axis([-3,3,-2,2,-0.7,1.5])[x1,y1]=meshgrid(-3:.2:3,-2:.2:2);z1=interp2(x,y,z,x1,y1);subplot(2,2,2);surf(x1,y1,z1),axis([-3,3,-2,2,-0.7,1.5])z2=interp2(x,y,z,x1,y1,'cubic');z3=interp2(x,y,z,x1,y1,'spline');subplot(2,2,3);surf(x1,y1,z2),axis([-3,3,-2,2,-0.7,1.5])subplot(2,2,4);surf(x1,y1,z3),axis([-3,3,-2,2,-0.7,1.5])画在同一幅图上作比较

前两个差值算法误差的比较z=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);surf(x1,y1,abs(z-z1)),axis([-3,3,-2,2,0,0.08])figure;surf(x1,y1,abs(z-z2)),axis([-3,3,-2,2,0,0.025])三个差值算法误差的比较z=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);subplot(3,1,1);surf(x1,y1,abs(z-z1)),axis([-3,3,-2,2,0,0.18])subplot(3,1,2);surf(x1,y1,abs(z-z2)),axis([-3,3,-2,2,0,0.18])subplot(3,1,3);surf(x1,y1,abs(z-z3)),axis([-3,3,-2,2,0,0.18])（2）二维一般分布数据的插值问题griddata函数的调用格式：z=griddata(x0,y0,z0,x,y,method)

method=‘v4’：插值算法，公认效果较好临近点插值：method=‘nearest’线性插值：method=‘linear’（缺省算法）三次样条插值：method=‘spline’立方插值：method=‘cubic’例9

在x为［-3,3]，y为［-2，2］矩形区域随机选择一组数据点，用’v4’与’cubic’插值法进行处理，并对误差进行比较。%产生数据点x=-3+6*rand(200,1);y=-2+4*rand(200,1);z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y);plot(x,y,'x')%样本点的二维分布

figure,plot3(x,y,z,'x'),axis([-3,3,-2,2,-0.7,1.5]),grid%cubic和v4算法[x1,y1]=meshgrid(-3:0.2:3,-2:0.2:2);z1=griddata(x,y,z,x1,y1,'cubic');subplot(2,2,1);surf(x1,y1,z1);axis([-3,3,-2,2,-0.7,1.5])z2=griddata(x,y,z,x1,y1,'v4');subplot(2,2,2);surf(x1,y1,z2);axis([-3,3,-2,2,-0.7,1.5])%误差分析z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);subplot(2,2,3);surf(x1,y1,abs(z0-z1)),axis([-3,3,-2,2,0,0.15])subplot(2,2,4);surf(x1,y1,abs(z0-z2));axis([-3,3,-2,2,0,0.15])4、高维插值问题三维插值interp3函数的调用格式：三维网格[x,y,z]=meshgrid(x1,y1,z1)griddata3()三维非网格数据插值n维插值interpn函数n维网格[x1,x2,…,xn]=ndgrid[v1,v2,…,vn]griddatan()n维非网格数据插值interp3()、interpn()调用格式同interp2()一致griddata3()、griddatan()调用格式同griddata()一致例10

通过函数生成一些网格型样本点，据此进行插值并给出插值误差。[x,y,z]=meshgrid(-1:0