范文一:数学建模插值与拟合
数据插值与拟合
插值与插值函数:已知由
的一批离散数据
(可能未知或非常复杂) 产生,且
个互异插值节点
,在插值区间内寻找一个相对简单的函
数
,使其满足下列插值条件:
再利用已求得的
计算任一非插值节点
的近似值
称为
,这就是插值。其中
被插函数。
称为插值函数,
最小二乘拟合: 已知一批离散的数据
互不相同,寻求一个拟合函数
,使
与
,
的误差平
方和在最小二乘意义下最小。在最小二乘意义下确定的
称为最小二乘拟合函数。 1)Lagrange 插值法
a .待定系数法: 假设插值多项式
,利用待定系数法即可求得满足插
值条件
数
。
的插值函数。关键在于确定待定系
b .利用基函数的构造方法 首先构造
的 次插值基函数
下的Lagrange 插值多项式:
个满足条件:
,再将其线性组合即可得如
其中
c .Lagrange 插值余项
注:上述两种构造方法所得的Lagrange 插值多项式是一样的,即满足插值条件
项式是唯一的。
2)分段线性插值
作分段线性插值的目的在于克服Lagrange 插值方法可能发生的不收敛性缺点。所谓分段线性插值就是利用每两个相邻插值节点作线性插值,即可得如下分段线性插值函数:
的Lagrange 插值多
其中
特点:插值函数序列
具有一致收敛性,克服了高次
Lagrange 插值方法的缺点,故可通过增加插值节点的方法提高其插值精度。但存在于节点处不光滑、插值精度低的缺点。 3)三次样条插值
三次样条插值的目的在于克服Lagrange 插值的不收敛性和提高分段线性插值函数在节点处的光滑性。所谓三次样条插值方法就是在满足下列条件: a
. b
.
在每个子区间
上是三次多项式的
三次样条函数中寻找满足如下插值条件:
以及形如
法。
特点:三次样条插值函数序列
一致收敛于被插函数,
等边界条件的插值函数
的方
因此可通过增加节点的方法提高插值的精度。
4)插值方法的Matlab 实现 一维数据插值
MATLAB 中用函数interp1来拟合一维数据,语法是 YI = INTERP1(X,Y,XI,方法)
其中(X , Y ) 是已给的数据点, XI 是插值点, 其中方法主要有
'linear' -线性插值,默认 'pchip' -逐段三次Hermite 插值 'spline' -逐段三次样条函数插值 其中最后一种插值的曲线比较平滑 例:
x=0:.12:1; x1=0:.02:1;
y=(x.^2-3*x+5).*exp(-5*x).*sin(x); plot(x,y,'o'); hold on; y1=interp1(x,y,x1,'spline'); plot(x1,y1,':')
如果要根据样本点求函数的定积分,而函数又是比较光滑的,则可以用样条函数进行插值后再积分,在MATLAB 中可以编写如下程序: function y=quadspln(x0,y0,a,b)
f=inline(‘interp1(x0,y0,x,’’spline ’’) ’, ’x ’, ’x0’, ’y0’); y=quadl(f,a,b,1e-8,[],x0,y0);
现求six(x)在区间[0,pi]上的定积分,只取5点 x0=[0,0.4,1,2,pi]; y0=sin(x0);
I=quadspln(x0,y0,0,pi)
结果得到的值为 2.01905, 精确值为2
二元函数插值:
MATLAB 中用函数interp2来拟合二维网格(X ,Y )上的数据Z ,语法是
YI = INTERP2(X,Y, Z,XI, YI,方法)
其中(X , Y ,Z ) 是已给的数据点, (XI ,YI )是插值点坐标, 其中方法主要有
'linear' -线性插值,默认 'pchip' -逐段三次Hermite 插值 'spline' -逐段三次样条函数插值 其中最后一种插值的曲面比较平滑
例:[x,y]=meshgrid(-3:.6:3,-2:.4:2); z=(x.^2-2*x).*exp(-x.^2-y.^2-x..*y);
[x1,y1]=meshgrid(-3:.2:3,-2:.2:2); %生成网格, x1和y1均为同样size 的矩阵
z1=interp2(x,y,z,x1,y1,’spline ’); %z1是矩阵,size 和x1,y1相同
surf(x1,y1,z1); axis([-3,3,-2,2,-0.7,1.5]);
-3
3
如果数据不是在网格上取的, 则可用函数griddata 来解决 语法是
YI = griddata(X,Y, Z,XI, YI,‘v4’)
其中(X , Y ,Z ) 是已给的数据点, (XI ,YI )是插值点坐标,
其中除了方法‘v4’外还有 'linear' -线性插值,默认 'cublc' -逐段三次Hermite 插值 'nearest' 其中‘v4’方法比较好
例x=-3+6*rand(200,1); %生成随机点的x 坐标向量x y=-2+4*rand(200,1); %生成随机点的y 坐标向量y z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); % 上述点的样本值向量z
[x1,y1]=meshgrid(-3:.2:3,-2:.2:2); %生成网格, x1和y1均为同样size 的矩阵 z1=griddata(x,y,z,x1,y1,’v4’); surf(x1,y1,z1); axis([-3,3,-2,2,-0.7,1.5]); 生成的图类似上图。
最小二乘拟合
假设有一组数据,xi, yi, 且已知它们满足某一函数y=f(a,x),其中a 是待定的参数,最小二乘拟合就是要确定这些参数,使得目标函数 sum(yi-f(a,xi))^2) 达到最小。在MATLAB 的最优化工具箱中有函数lsqcurvefit, 语法如下
[a,J]=lsqcurvefit(fun,a0,x,y,Lb,Ub,options), 其中fun 是内联函数,也可以是m 函数,但要写为@fun的形式, a0是参数的初值,x,y 为 数据向量。J=min( sum(yi-f(a,xi))^2)) Lb 和Ub 是a 的下上限 例 x=0:.1:10;
y=0.12*exp(-0.213*x)+0.54*exp(-0.17*x).*sin(1.23*x); f=inline('a(1)*exp(-a(2)*x)+a(3)*exp(-a(4)*x).*sin(a(5)*x)','a','x');
[xx,res]=lsqcurvefit(f,[1,1,1,1,1],x,y),
运行得a 向量为 0.1197 0.2125 0.5404 0.1702 1.2300
最小二乘为 7.1637e-007
如果要提高精度, 可把上面程序的最后一句改为以下语句。首先自定义选项 Ff=optimset;
ff.TolFun=1e-20; ff.TolX=1e-15,
[xx,res]=lsqcurvefit(f,[1,1,1,1,1],x,y,[],[],ff);
题目:拟合问题
在农业生产试验研究中,对某地区土豆的产量与化肥的关系做了一实验,得到了氮肥、磷肥的施肥量与土豆产量的对应关系如下表:
根据上表数据分别给出土豆产量与氮、磷肥的关系式。
建模分析:
使用 Matlab 语言首先画出土豆产量与氮施肥量的散点图,从图可看出土豆产量与氮肥量的关系是二次函数关系,因此可选取拟合函数为:
,
其中 x 和 y 分别为氮肥量和土豆产量,a 、b 和 c 为待定系数。再画出磷肥量与土豆产量的散点图,从图可看出从 0 到 98、从 98 到 342 之间分别呈明显的线性关系。由此可选取所求拟合函数为一分段的线性函数,换言之,用前 5 点作一线性拟合函数,再用后 6 个点也作一线性拟合函数,最后用两个线性函数求出其分界点即可得分段线性函数。 数学模型:
对氮肥的拟合函数为:
对磷肥的拟合函数为:
x1=[0,34,67,101,135,202,259,336,404,471];
y1=[15.18,21.36,25.72,32.29,34.03,39.45,43.15,43.46,40.83,30.75]; plot(x1,y1,'r+') aa=polyfit(x1,y1,2) xx=0:471;
yy=aa(1)*xx.*xx+aa(2)*xx+aa(3); plot(xx,yy,x1,y1,'r+') %Next problem
x2=[0,24,49,73,98,147,196,245,294,342];
y2=[33.46,32.47,36.06,37.96,41.04,40.09,41.26,42.17,40.36,42.73]; plot(x2,y2,'r+')
a1=polyfit(x2(1:5),y2(1:5),1) a2=polyfit(x2(5:10),y2(5:10),1) x0=(a2(2)-a1(2))/(a1(1)-a2(1)) xx1=0:x0; yy1=a1(1)*xx1+a1(2); xx2=x0:342; yy2=a2(1)*xx2+a2(2); plot(x2,y2,'r+',xx1,yy1,xx2,yy2) 程序运行结果
11
范文二:数学建模案例分析-- 插值与拟合方法建模1数据插值方法及应用
第十章 插值与拟合方法建模
在生产实际中~常常要处理由实验或测量所得到的一批离散数据~插值与拟合方法就是要通过这些数据去确定某一类已经函数的参数~或寻求某个近似函数使之与已知数据有较高的拟合精度。插值与拟合的方法很多~这里主要介绍线性插值方法、多项式插值方法和样条插值方法~以及最小二乘拟合方法在实际问题中的应用。相应的理论和算法是数值分析的内容~这里不作详细介绍~请参阅有关的书籍。
?1 数据插值方法及应用
在生产实践和科学研究中,常常有这样的问题:由实验或测量得到变量间的一批离散样点,要求由此建立变量之间的函数关系或得到样点之外的数据。与此有关的一类问题是当原始数据
精度较高,要求确定一个初等函数(一般用多项式或分段(x,y),(x,y),?,(x,y)y,P(x)0011nn
多项式函数)通过已知各数据点(节点),即,或要求得函数在另外一y,P(x),i,0,1,?,nii
些点(插值点)处的数值,这便是插值问题。
1、分段线性插值
这是最通俗的一种方法,直观上就是将各数据点用折线连接起来。如果
a,x,x,?,x,b01n
那么分段线性插值公式为
x,xx,xii,1P(x),y,y,x,x,x,i,1,2,?,n i,1ii,1ix,xx,xi,1iii,1
可以证明,当分点足够细时,分段线性插值是收敛的。其缺点是不能形成一条光滑曲线。 例1、已知欧洲一个国家的地图,为了算出它的国土面积,对地图作了如下测量:以由西向东方向为x轴,由南向北方向为y轴,选择方便的原点,并将从最西边界点到最东边界点在x轴上的区间适当的分为若干段,在每个分点的y方向测出南边界点和北边界点的y坐标y1和y2,这样就得到下表的数据(单位:mm)。
x 7.0 10.5 13.0 17.5 34.0 40.5 44.5 48.0 56.0 y1 44 45 47 50 50 38 30 30 34 y2 44 59 70 72 93 100 110 110 110 x 61.0 68.5 76.5 80.5 91.0 96.0 101.0 104.0 106.5 y1 36 34 41 45 46 43 37 33 28 y2 117 118 116 118 118 121 124 121 121 x 111.5 118.0 123.5 136.5 142.0 146.0 150.0 157.0 158.0 y1 32 65 55 54 52 50 66 66 68 y2 121 122 116 83 81 82 86 85 68
根据地图的比例,18 mm相当于40 km。
根据测量数据,利用MATLAB软件对上下边界进行线性多项式插值,分别求出上边界函数,f(x)2下边界函数,利用求平面图形面积的数值积分方法—将该面积近似分成若干个小长方形,f(x)1
分别求出这些长方形的面积后相加即为该面积的近似解。
n
S,lim[f(,),f(,)],x,2i1ii,,n,i1
式中,。 ,,[x,x]ii,1i
这里线性插值和面积计算源程序如下:
clear all
x=[7.0 10.5 13.0 17.5 34.0 40.5 44.5 48.0 56.0 61.0 68.5 76.5 80.5 91.0 96.0 101.0 104.0 106.5 111.5
118.0 123.5 136.5 142.0 146.0 150.0 157.0 158.0]; y1=[44 45 47 50 50 38 30 30 34 36 34 41 45 46 43 37 33 28 32 65 55 54 52 50 66 66 68];
y2=[44 59 70 72 93 100 110 110 110 117 118 116 118 118 121 124 121 121 121 122 116 83 81 82 86 85
68];
newx=7:0.1:158;
newy1=interp1(x,y1,newx,’linear’);
newy2=interp1(x,y2,newx,’linear’);
Area=sum(newy2- newy1)*0.1/18^2*1600 最后计算的面积约为42414平方公里。
2、多项式插值
设有次多项式 m
mm,1 P(x),ax,ax,?,ax,a01m,1m
n,1通过所有个点(x,y),(x,y),?,(x,y),那么就有 0011nn
mm,1ax,ax,?,ax,a,y,i,0,1,?,n 0i1im,1imi
可以证明当x,x,?,x且时,这样的多项式存在且唯一。若要求得到函数表达式,m,n01n
可直接解上面方程组。若只要求得函数在插值点处数值,可用下列Lagrange插值公式
nnx,**x,y()() ,,nix,x,0,0,,ijjiij
m,7多项式插值光滑但不具有收敛性,一般不宜采用高次多项式(如)插值。 例2、在万能拉拨机中有一个园柱形凸轮,其底园半径R=300mm,凸轮的上端面不在同一平面上,
i而要根据动杆位移变化的需要进行设计制造。按设计要求,将底园周18等分,旋转一周。第个
y(i,0,1,2,?,18)分点对应柱高,数据见下表。为了数控加工,需要计算出园周上任一点的柱i
高。
凸轮高度的数据,单位:mm,
i分点 0和18 1 2 3 4 5
柱高 502.8 525.0 514.3 451.0 326.5 188.6
i分点 6 7 8 9 10 11
柱高 92.2 59.6 62.2 102.7 147.1 191.6
i分点 12 13 14 15 16 17
柱高 236.0 280.5 324.9 369.4 413.8 458.3
我们将园周展开,借助MATLAB软件画出对应的柱高曲线散点图(左下图)。 clear;close;
x=linspace(0,2*pi*300,19);
y=[502.8 ,525.0,514.3,451.0,326.5,188.6,92.2,59.6,62.2,102.7,147.1,191.6,236.0,280.5,324.9,369.4,413.
8,458.3,502.8];
plot(x,y,’o’);axis([0,2000,0,550]);
可见,可以用三次多项式插值,下面给出借助MATLAB软件画出的柱高插值曲线图(右上图)。 xi=0:2*pi*300;
yi=interp1(x,y,xi,’cubic’);
plot(xi,yi);
3、样条插值
这是最常用的插值方法。数学上所说的样条,实质上是指分段多项式的光滑连接。设有
a,x,x,?,x,b 01n
称分段函数为k次样条函数,若它满足 S(x)
k(1) 在每个小区间上是次数不超过次的多项式; S(x)
k,1(2) 在上具有直到阶的连续导数。 S(x)[a,b]
用样条函数作出的插值称为样条插值。工程上广泛采用三次样条插值。 例3、某居民区的自来水是由一个园柱形的水塔提供。水塔高12.2米,直径17.4米。水塔由水泵根据塔中水位高低自动加水,一般每天水泵工作两次。按照设计,当水塔内的水位降至约8.2米时,水泵自动启动加水;当水位升至约10.8米时,水泵停止工作。现在需要了解该居民区用水规律,这可以通过用水率(单位时间的用水量)来反映。通过间隔一段时间测量水塔中的水位来估算用水率。下表是某一天的测量记录数据,测量了28个时刻(单位:小时)的水位(单位:米),但由于其中有3个时刻正遇到水泵在向水塔供水,而无水位记录(表中用符号//表示)。
时刻 0 0.921 1.843 2.949 3.871 4.978 5.900
水位 9.677 9.479 9.308 9.125 8.982 8.814 8.686
时刻 7.006 7.982 8.967 9.981 10.925 10.954 12.032
水位 8.525 8.388 8.220 // // 10.820 10.500
时刻 12.954 13.875 14.982 15.903 16.826 17.931 19.037
水位 10.210 9.936 9.653 9.409 9.180 8.921 8.662
时刻 19.959 20.839 22.015 22.958 23.880 24.986 25.908
水位 8.433 8.220 // 10.820 10.591 10.354 10.180
,2h先通过体积公式,利用上表中的水位高,得到不同时刻水塔中水的体积。为提v,dhtvii4
2高精度,采用二阶差商来估算时刻的水流速度,即。 tf(t),,,viii
具体地,因为所有数据被水泵两次工作分割成三组数据,对每组数据的中间数据采用中心差商,前后两个数据不能够采用中心差商,改用向前或向后差商。
,v,8v,8v,v2i,2i,1i,1i,2中心差商公式 ,v, i12(t,t)i,1i
,v,4v,3v2i,2i,1i,v,向前差商公式 i2(t,t)i,1i
3v,4v,v2ii,1i,2,v,向后差商公式 i2(t,t)ii,1
估算出水塔中水的流速(单位:立方米/小时)见下表。
时刻 0 0.921 1.843 2.949 3.871 4.978 5.900
流速 54.516 42.320 38.085 41.679 33.297 37.814 30.748
时刻 7.006 7.982 8.967 9.981 10.925 10.954 12.032
流速 38.455 32.122 41.718 // // 73.686 76.434
时刻 12.954 13.875 14.982 15.903 16.826 17.931 19.037
流速 71.686 60.190 68.333 59.217 52.011 56.626 63.023
时刻 19.959 20.839 22.015 22.958 23.880 24.986 25.908
流速 54.859 55.439 // 57.602 57.766 51.891 36.464
先用MATLAB画出水流速散点图。
t=[0 0.921 1.843 2.949 3.871 4.978 5.9 7.006 7.982 8.967 10.954 12.032 12.954 13.875 14.982 15.903
16.826 17.931 19.037 19.959 20.839 22.958 23.88 24.986 25.908];
r=[54.516 42.320 38.085 41.679 33.297 37.814 30.748 38.455 32.122 41.718 73.686 76.434 71.686
60.19 68.333 59.217 52.011 56.626 63.023 54.859 55.439 57.602 57.766 51.891 36.464];
plot(t,r,’b+’); % (t,r)表示时间和流速
title(‘流速散点图’);xlabel(’时间(小时)’); ylabel(‘流速(立方米/小时)’)
使用MATLAB软件中的三次样条插值命令得到用水率函数f(t)如下图所示。 x0=t;y0=r;
[l,n]=size (x0); dl=x0(n)-x0(1);
x=x0(1):1/3600:x0(n); %被插值点
ys=interp1 (x0,y0,x,’spline’); %样条插值输出
plot (x,ys);
title(‘样条插值下的流速图’);xlabel(’时间(小时)’); ylabel(‘流速(立方米/小时)’)
下面是余秋雨经典励志语录,欢迎阅读。
不需要的朋友可以编辑删除~~
关于年龄
1.一个横贯终生的品德基本上都是在青年时代形成的,可惜在那个至关重要的时代,青年人受到的正面的鼓动永远是为成功而搏斗,而一般所谓的成功总是带有排他性、自私性的印记。结果,脸颊上还没有皱纹的他们,却在品德上挖下了一个个看不见的黑洞。
2.我不赞成太多地歌颂青年,而坚持认为那是一个充满陷阱的年代。陷阱一生都会遇到,但青年时代的陷阱最多、最大、最险。
3.历史上也有一些深刻的哲人,以歌颂青年来弘扬社会的生命力。但这里显然横亘着一种二律背反:越是坚固的对象越需要鼓动青年去对付,但他们恰恰因为年轻,无法与真正的坚持相斡旋。
4.青年时代的正常状态是什么,我想一切还是从真诚的谦虚开始。青年人应该懂得,在我们出生之前,这个世界已经精精彩彩、复复杂杂地存在过无数年,我们什么也不懂,能够站正脚下的一角建设一点什么,已是万幸。
5.中年是对青年的延伸,又是对青年的告别。这种告别不仅仅是一系列观念的变异,而是一个终于自立的成熟者对于能够随心所欲处置各种问题的自信。
6.中年人的当家体验是最后一次精神断奶。你突然感觉到终于摆脱了父母、兄长、老师的某种依赖,而这种依赖在青年时代总是依稀犹在的;对于领导和组织,似乎更贴近了,却又显示出自己的独立存在,你成了社会结构网络中不可缺少的一个点;因此你在热闹中品尝了有生以来真正的孤立无援,空前的脆弱和空前的强大集于一身。
7.中年人一旦有了当家体验,就会明白教科书式的人生教条十分可笑。当家管着这么一个大摊子,每个角落每时每刻都在涌现着新问题,除了敏锐而又细致地体察实际情况,实事求是地解开每一个症结,简直没有高谈阔论、把玩概念的余地。这时人生变得很空灵,除了隐隐然几条人生大原则,再也记不得更多的条令。
8.中年人的坚守,已从观点上升到人格,而人格难以言表,他们变得似乎已经没有顶在脑门上的观点。他们知道,只要坚守着自身的人格原则,很多看似对立的观点都可相容相依,一一点化成合理的存在。于是,在中年人眼前,大批的对峙消解了,早年的对手找不到了,昨天的敌人也没有太多仇恨了,更多的是把老老少少各色人等照顾在自己身边。请不要小看这“照顾”二字,中年人的魅力至少有一半与此相关。
9.中年人最可怕的是失去方寸。这比青年人和老年人的失态有更大的危害。中年人失去方寸的主要特征是忘记自己的年龄。一会儿要别人像对待青年那样关爱自己,一会儿又要别人像对待老人那样尊敬自己,他永远活在中年之外的两端,偏偏不肯在自己的年龄里落脚。
10、某个时期,某个社会,即使所有的青年人和老年人都中魔一般荒唐了,只要中年人不荒唐,事情就坏不到哪里去。最怕的是中年人的荒唐,而中年人最大的荒唐,就是忘记了自己是中年。
11、中年太实际、太繁忙,在整体上算不得诗,想来难理解;青年时代常常被诗化,但青年时代的诗太多激情而缺少意境,按我的标准,缺少意境就算不得好诗。
12、一般情况下,老年岁月总是比较悠闲,总是能够没有功利而重新面对自然,总是漫步在回忆的原野,而这一切,都是诗和文学的特质所在。老年人可能不会写诗或已经不再写诗,但他们却以诗的方式生存着。看街市忙碌,看后辈来去,看庭花凋零,看春草又绿,而思绪则时断时续,时喜时悲,时真时幻。 13、老人的年龄也有积极的缓释功能,为中青年的社会减轻负担。不负责任的中青年用不正当的宠溺败坏了老人的年龄,但老人中毕竟还有冷静的智者,默默固守着年岁给予的淡然的尊严。
14、只有到了老年,沉重的人生使命已经卸除,生活的甘苦也已了然,万丈红尘已移到远处,宁静下来了的周际环境和逐渐放慢了的生命节奏构成了一种总结性、归纳性的轻微和声,诗的意境出现了。
15、中青年的世界再强悍,也经常需要一些苍老的手来救助。平时不容易见到,一旦有事则及时伸出,救助过后又立即消失,神龙见首不见尾。这是一种早已退出社会主体的隐性文化和柔性文化,隐柔中沉积着岁月的硬度,能使后人一时启悟,如与天人对晤。老年的魅力,理应在这样的高位上偶尔显露。不要驱使,不要强求,不要哄抬,只让它们成为人生的写意笔墨,似淡似浓,似有似无。
关于人生
1.我们对这个世界,知道得还实在太少。无数的未知包围着我们,才使人生保留迸发的乐趣。当哪一天,世界上的一切都能明确解释了,这个世界也就变得十分无聊。人生,就会成为一种简单的轨迹,一种沉闷的重复。
2.人有多种活法,活着的文明等级也不相同,住在五层楼上的人完全不必去批评三层楼的低下,何况你是否在五层楼还缺少科学论证。
3.人生的道路也就是从出生地出发,越走越远。一出生便是自己,由此开始的人生就是要让自己与种种异己的一切打交道。打交道的结果可能丧失自己,也可能在一个更高的层面上把自己找回。
4.不管你今后如何重要,总会有一天从热闹中逃亡,孤舟单骑,只想与高山流水对晤。走得远了,也许会遇到一个人,像樵夫,像路人,出现在你与高山流水之间,短短几句话,
使你大惊失色,引为终生莫逆。但是,天道容不下如此至善至美,你注定会失去他,同时也就失去了你的大半生命。
5.人生的过程虽然会受到社会和时代的很大影响,但贯穿首尾的基本线索总离不开自己的个体生命。个体生命的完整性、连贯性会构成一种巨大的力量,使人生的任何一个小点都指点着整体价值。
6.如果有一天,我们突然发现,投身再大的事业也不如把自己的人生当做一个事业,聆听再好的故事也不如把自己的人生当做一个故事,我们一定会动手动笔,做一点有意思的事情。
7.杰出之所以杰出,是因为罕见,我们把自己连接于罕见,岂不冒险?既然大家都很普通,那么就不要鄙视世俗岁月、庸常岁序。不孤注一掷,不赌咒发誓,不祈求奇迹,不想入非非,只是平缓而负责地一天天走下去,走在记忆和向往的双向路途上,这样,平常中也就出现了滋味,出现了境界。
8.就人生而言,应平衡于山、水之间。水边给人喜悦,山地给人安慰。水边让我们感知世界无常,山地让我们领悟天地恒昌。水边让我们享受脱离长辈怀抱的远行刺激,山地让我们体验回归祖先居所的悠悠厚味。
9.第一根白发人人都会遇到,谁也无法讳避,因此这个悲剧似小实大,简直是天网恢恢,疏而不漏,而决斗、毒药和暗杀只是偶发性事件,这种偶发性事件能快速置人于死地,但第一根白发却把生命的起点和终点连成了一条绵长的逻辑线,人生的任何一段都与它相连。
10、谁也不要躲避和掩盖一些最质朴、最自然的人生课题如年龄问题。再高的职位,再多的财富,再大灾难,比之于韶华流逝、岁月沧桑、长幼对视、生死交错,都成了皮相。北雁长鸣,年迈的帝王和年迈的乞丐一起都听到了;寒山扫墓,长辈的泪滴和晚辈的泪滴却有不同的重量。
11、人格尊严的表现不仅仅是强硬。强硬只是人格的外层警卫。到了内层,人格的天地是清风明月,柔枝涟漪,细步款款,浅笑连连。
12、黄山谷说过:“人胸中久不用古今浇灌,则尘俗生其间,照镜觉面目可憎,对人亦语言无味。”这就是平庸的写照。如此好事,如果等到成年后再来匆匆弥补就有点可惜了,最好在青年时就进入。早一天,就多一份人生的精彩;迟一天,就多一天平庸的困扰。
13、再高的职位,再多的财富,再大灾难,比之于韶华流逝、岁月沧桑、长幼对视、生死交错,都成了皮相。北雁长鸣,年迈的帝王和年迈的乞丐一起都听到了;寒山扫墓,长辈的泪滴和晚辈的泪滴却有不同的重量。
14、人生不要光做加法。在人际交往上,经常减肥、排毒,才会轻轻松松地走以后的路。
15、几乎每一个改革探索者都遇到过嫉妒的侵扰,更不要说其中的成功者了。人们很容易对高出自己视线的一切存在投去不信任,在别人快速成功的背后寻找投机取巧的秘密。
关于文化
1.真正的文化精英是存在的,而且对国家社会非常重要。但是这些年来,由于伪精英的架势实在是太让人恶心了,结果连真的精英的名声也败坏了。真精英总是着眼于责任,伪精英总是忙着装扮;真精英总是努力地与民众沟通,伪精英总是努力地与民众划分,这就是最根本的区别。
2.凡是文化程度不高的群落,总是会对自己不懂的文化话语心存敬畏,正是这种敬畏心理被一些投机文人利用了。
3.在文化上,无效必然导致无聊,无聊又必然引来无耻。但是,即使到了这种“三无”的低谷,也不必过于沮丧。因为只有低谷,才能构成对新高峰的向往。
4.当今天下百业,文化最大。当今天下百行,文化届最小。那么,岂能再让一个日渐干涸的小池塘,担任江河湖海的形象代表?
5.古代绘画中无论是萧瑟的荒江、丛山中的苦旅,还是春光中的飞鸟、危崖上的雏鹰,只要是传世佳品,都会包藏着深厚的人生意识。贝多芬的交响曲,都是人生交响曲。
6.善良,这是一个最单纯的词汇,又是一个最复杂的词汇。它浅显到人人都能领会,又深奥到无人能够定义。它与人终生相伴,但人们却很少琢磨它、追问它。
7.社会理性使命已悄悄抽绎,秀丽山水间散落着才子、隐士,埋藏着身前的孤傲和身后的空名。天大的才华和郁愤,最后都化作供后人游玩的景点。
8.阅读的最大理由是想摆脱平庸,早一天就多一份人生的精彩;迟一天就多一天平庸的困扰。
9.为什么那么多中国民众突然对韩国的电视剧,对超女表现出那么单纯的投入,很重要的原因是,韩国艺术家不知道中国评论家,而超女根本不在乎评论家的存在。
10、一切美丽都是和谐的,因此总是浑然天成,典雅含蓄。反之,一切丑陋都是狞厉的,因此总是耀武扬威,嚣张霸道。如果没有审美公德的佑护,美永远战胜不了丑。
11、什么季节观什么景,什么时令赏什么花,这才完整和自然。如果故意地大颠大倒,就会把两头的况味都损害了。“暖冬”和“寒春”都不是正常的天象。
12、文明的人类总是热衷于考古,就是想把压缩在泥土里的历史扒剔出来,舒展开来,窥探自己先辈的种种真相。那么,考古也就是回乡,也就是探家。探视地面上的家乡往往会有岁月的唏嘘、难言的失落,使无数游子欲往而退;探视地底下的家乡就没有那么多心理障碍了,整个儿洋溢着历史的诗情、想像的愉悦。
13、我们的历史太长、权谋太深、兵法太多、黑箱太大、内幕太厚、口舌太贪、眼光太杂、预计太险,因此,对一切都“构思过度”。
14、中华文化的三大优点:一、不喜远征。中国人不会举一国之力去攻打远方之国。二、不喜极端。儒家讲究“中庸之道”,会努力寻找一个中间点,规避极端三、不喜无序。中国一直处于集权统治的状态中,习惯所有的事务都在管理之中,中国失控的时候是很少见的。
关于爱情
很多女孩子觉得责任感不太重要,男人没有责任感反而给了女方一种权利。其实对男人来说,还有什么比没有责任感可怕地呢?与没有责任感的男人谈恋爱,就像与朝雾和晚霞厮磨,再美好也没有着落。
爱情非常珍贵,不仅值得用斗争来保卫,而且即使付出生命的代价也值得。
其实,未经艰苦寻找的草率结合,对她也是不尊重。她和你一样,都有寻求深刻爱情的权利。
每一男女都处在自转之中,当一个男人最散发魅力的一面转向了一位女人,而这女人最美好的一面也刚好朝向了这个男人,那么爱情就挡也挡不住了。当然不是每个人都如此幸运,自转的方向和速度,相对于那个有可能出现或已经错过的异性,总要有偏差,所以老有人找不到自己的爱情。
2、能够慢慢培养的不是爱情,而是习惯。能够随着时间得到的,不是感情而是感动。所以爱是一瞬间的礼物,有就有,没有就没有。但反过来说,爱和婚姻实际并不是一回事情,并不是所有的爱情都要结婚的,也不是所有婚姻都有爱情的。
6、爱情里,总有一个主角和一个配角,累的永远是主角,伤的永远是配角;有时,爱也是种伤害:残忍的人,选择伤害别人,善良的人,选择伤害自己;人生就是一种承受,需要学会支撑。支撑事业,支撑家庭,甚至支撑起整个社会,有支撑就一定会有承受,支撑起多少重量,就要承受多大压力。
7、假如你想要一件东西,就放它走。它若能回来找你,就永远属于你;它若不回来,那根本就不是你的。爱情也是如此。
8、为什么把择定终身的职责,交付给半懂不懂的年岁;为什么把成熟的眼光,延误地出现在早已收获过的荒原?
9、说了那么多旳——“如若你不在,我等待你归来。”也比不过你一句——“我不会等,我去找你!”
关于友情
1.常听人说,人世间最纯净的友情只存在于孩童时代。这是一句极其悲凉的话,居然有那么多人赞成,人生之孤独和艰难,可想而知。我并不赞成这句话。孩童时代的友情只是愉快的嘻戏,成年人靠着回忆追加给它的东西很不真实。友情的真正意义产生于成年之后,它不可能在尚未获得意义之时便抵达最佳状态。
2.很多人都是在某次友情感受的突变中,猛然发现自己长大的。仿佛是哪一天的中午或傍晚,一位要好同学遇到的困难使你感到了一种不可推卸的责任,你放慢脚步忧思起来,开始懂得人生的重量。就在这一刻,你突然长大。
3.在人生的诸多荒诞中,首当其冲的便是友情的错位。友情的错位,来源于我们自身的混乱。
4.置身于同一个职业难道是友情的基础?当然不是。如果偶尔有之,也不能本末倒置。情感岂能依附于事功,友谊岂能从属于谋生,朋友岂能局限于同僚。
5.在家靠父母,出外靠朋友。这种说法既表明了朋友的重要,又表明了朋友的价值在于被依靠。但是,没有可靠的实用价值能不能成为朋友?一切帮助过你的人是不是都能算作朋友?
6.患难见知己,烈火炼真金。这又对友情提出了一种要求,盼望它在危难之际及时出现。能够出现当然很好,但友情不是应急的储备,朋友更不应该被故意地考验。
7.真正的友情不依靠什么。不依靠事业、祸福和身份,不依靠经历、方位和处境,它在本性上拒绝功利,拒绝归属,拒绝契约,它是独立人格之间的互相呼应和确认。它使人们独而不孤,互相解读自己存在的意义。因此所谓朋友也只不过是互相使对方活得更加自在的那些人。
8.真正的友情都应该具有“无所求” 的性质,一旦有所求,“求”也就成了目的,友情却转化为一种外在的装点。我认为,世间的友情至少有一半是被有所求败坏的,即便所求的内容乍一看并不是坏东西;让友情分担忧愁,让友情推进工作??,友情成了忙忙碌碌的工具,那它自身又是什么呢?应该为友情卸除重担,也让朋友们轻松起来。朋友就是朋友,除此之外,无所求。
9.无所求的朋友最难得,不妨闭眼一试,把有所求的朋友一一删去,最后还剩几个?
10.真正的友情因为不企求什么不依靠什么,总是既纯净又脆弱。世间的一切孤独者也都遭遇过友情,只是不知鉴别和维护,一一破碎了。
11.“君子之交谈如水”,这种高明的说法包藏着一种机智的无奈,可惜后来一直被并无机智、只剩无奈的人群所套用。怕一切许诺无法兑现,于是不作许诺;怕一切欢晤无法延续,于是不作欢晤,只把微笑点头维系于影影绰绰之间。有人还曾经借用神秘的东方美学来支持这种态度:只可意会,不可言传;不着一字,尽得风流;羚羊挂角,无迹可寻??这样一来,友情也就成了一种水墨写意,若有若无。但是,事情到了这个地步,友情和相识还有什么区别?
12.强者捆扎友情,雅者淡化友情,俗者粘贴友情,都是为了防范友情的破碎,但看来看去,没有一个是好办法。原因可能在于,这些办法都过分依赖技术性手段,而技术性手段一旦进入感情领域,总没有好结果。
13.万不能把防范友情的破碎当成一个目的。该破碎的让它破碎,毫不足惜;虽然没有破碎却发现与自己生命的高贵内质有严重羝牾,也要做破碎化处理。罗丹说,什么是雕塑?那就是在石料上去掉那些不要的东西。我们自身的雕塑,也要用力凿掉那些异己的、却以朋友名义贴附着的杂质。不凿掉,就没有一个像模像样的自己。
14.该破碎的友情常被我们捆扎、粘合着,而不该破碎的友情却又常常被我们捏碎了。两种情况都是悲剧,但不该破碎的友情是那么珍贵,它居然被我们亲手捏碎,这对人类良知的打击几乎是致命的。
15.其实,世上哪有两片完全相同的树叶,即便这两片树叶贴得很紧?本有差异却没有差异准备,都把差异当作了背叛,夸张其词地要求对方纠正。这是一种双方的委屈,友情的回忆又使这种委屈增加了重量。负荷着这样的重量不可能再来纠正自己,双方都怒气冲天地走上了不归路。凡是重友情、讲正气的人都会产生这种怒气,而只有小人才是不会愤怒的一群,因此正人君子们一旦落入这种心理陷阱往往很难跳得出来。高贵的灵魂吞咽着说不出口的细小原因在陷阱里挣扎。
16.友情好像是一台魔力无边的红外线探测仪,能把一切隐藏的角落照个明明白白。不明不白也不要紧,理解就是一切,朋友总能理解,不理解还算朋友?但是,当误会无可避免地终于产生时,原先的不明不白全都成了疑点,这对被疑的一方而言无异是冤案加身;申诉无门,他的表现一定异常,异常的表现只能引起更大的怀疑,互相的友情立即变得难于收拾。
17.友情本是超越障碍的翅膀,但它自身也会背负障碍的沉重,因此,它在轻松人类的时候也在轻松自己,净化人类的时候也在净化自己。其结果应该是两相完满:当人类在最深刻地享受友情时,友情本身也获得最充分的实现。
18.现在,即便我们拥有不少友情,它也还是残缺的,原因在于我们自身还残缺。世界理应给我们更多的爱,我们理应给世界更多的爱,这在青年时代是一种小心翼翼的企盼,到了生命的秋季,仍然是一种小心翼翼的企盼。但是,秋季毕竟是秋季,生命已承受霜降,企盼已洒上寒露,友情的渴望灿如枫叶,却也已开始飘落。
范文三:数学建模插值及拟合详解
插值和拟合
实验目的:了解数值分析建模的方法,掌握用Matlab 进行曲线拟合的方法,理解用插值法建模的思想,运用Matlab 一些命令及编程实现插值建模。
实验要求:理解曲线拟合和插值方法的思想,熟悉Matlab 相关的命令,完成相应的练习,并将操作过程、程序及结果记录下来。 实验内容: 一、插值
1.插值的基本思想
·已知有n +1个节点(xj,yj ),j = 0,1,…, n,其中xj 互不相同,节点(xj, yj)可看成由某个函数 y= f(x )产生;
·构造一个相对简单的函数 y=P(x);
·使P 通过全部节点,即 P (xk) = yk,k=0,1,…, n ; ·用P (x)作为函数f ( x )的近似。 2.用MA TLAB 作一维插值计算 yi=interp1(x,y ,xi ,'method')
注:yi —xi 处的插值结果;x ,y —插值节点;xi —被插值点;method —插值方法(‘nearest’ :最邻近插值;‘linear’ :线性插值;‘spline’ :三次样条插值;‘cubic’ :立方插值;缺省时:线性插值)。 注意:所有的插值方法都要求x 是单调的,并且xi 不能够超过x 的范围。 练习1:机床加工问题
每一刀只能沿x 方向和y 方向走非常小的一步。 表3-1给出了下轮廓线上的部分数据
但工艺要求铣床沿x 方向每次只能移动0.1单位. 这时需求出当x 坐标每改变0.1单位时的y 坐标。 试完成加工所需的数据,画出曲线. 步骤1:用x0,y0两向量表示插值节点;
步骤2:被插值点x=0:0.1:15; y=y=interp1(x0,y0,x,'spline'); 步骤3:plot(x0,y0,'k+',x,y,'r')
grid on
答:x0=[0 3 5 7 9 11 12 13 14 15 ];
y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6 ]; x=0:0.1:15;
y=interp1(x0,y0,x,'spline'); plot(x0,y0,'k+',x,y,'r') grid on
1
3.用MA TLAB 作网格节点数据的插值(二维) z=interp2(x0,y0,z0,x,y,’method’) 注:z —被插点值的函数值;x0,y0,z0—插值节点;x ,y —被插值点;method —插值方法(‘nearest’ :最邻近插值;‘linear’ :双线性插值; ‘cubic’ :双三次插值;缺省时:双线性插值)。 注意:要求x0,y0单调;x ,y 可取为矩阵,或x 取行向量,y 取为列向量,x,y 的值分别不能超出x0,y0的范围。
4.用MA TLAB 作散点数据的插值计算 cz =griddata(
x ,y ,
z ,cx ,cy ,
‘method’)
注:cz —被插点值的函数值;x ,y ,z —插值节点;cx ,cy —被插值点;method —插值方法(‘nearest’ :最邻近插值;‘linear’ :双线性插值; ‘cubic’ :双三次插值;'v4‘:Matlab 提供的插值方法;缺省时:双线性插值)。
练习2:航行区域的警示线
某海域上频繁地有各种吨位的船只经过。 为保证船只的航行安全,有关机构在低潮时对水深进行了测量,下表是他们提供的测量数据: 水道水深的测量数据
x 129.0 140.0 103.5 88.0 185.5 195.0 105.5 y 7.5 141.5 23.0 147.0 22.5 137.5 85.5 z 4 8 6 8 6 8 8
x 157.5 107.5 77.0 81.0 162.0 162.0 117.5 y -6.5 -81.0 3.0 56.5 -66.5 84.0 -33.5 z 9 9 8 8 9 4 9
2
其中(x, y)为测量点,z 为(x, y)处的水深(英尺) ,水深z 是区域坐标(x, y)的函数z= z (x, y ),
船的吨位可以用其吃水深度来反映,分为 4英尺、4.5英尺、5英尺和 5.5英尺 4 档。
航运部门要在矩形海域(75,200)×(-50,150)上为不同吨位的航船设置警示标记。
请根据测量的数据描述该海域的地貌,并绘制不同吨位的警示线,供航运部门使用。 x=[129 140 103.5 88 185.5 195 105.5 157.5 107.5 77 81 162 162 117.5]; y=[7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5]; z=[-4 -8 -6 -8 -6 -8 -8 -9 -9 -8 -8 -9 -4 -9];
cx=75:0.5:200; cy=-70:0.5:150;
cz=griddata(x,y,z,cx,cy','cubic');
meshz(cx,cy,cz),rotate3d
xlabel('X'),ylabel('Y'),zlabel('Z') %pause
figure(2),contour(cx,cy,cz,[-5 -5]);grid on, hold on plot(x,y,'+')
xlabel('X'),ylabel('Y')
Z
200
Y
50
X
200
3
140120
100806040200-20-40-60
80
100
120
140X
160
180
200
Y
练习3:估计水塔的水流量
问题描述见教材P.91—93,请绘出三次样条插值曲线,并计算一天的总的用水量。 解:t0=[0.46,1.38,2.4,3.41,4.43,5.44,6.45,7.47,8.45,11.49,12.49,13.42,14.43,15.44,16.37,17.38,18.49,19.50,20.40,24.43,25.32];
v0=[11.2,9.7,8.6,8.1,9.3,7.2,7.9,7.4,8.4,15.6,16.4,15.5,13.4,13.8,12.9,12.2,12.2,12.9,12.6,11.2,3.5]; t=0:0.1:26; y=interp1(t0,v0,t,'spline'); plot(t0,v0,'k+',t,y, 'r') grid on
4
20
15
10
5
-5
-10
051015202530
二、曲线拟合
已知一组(二维)数据,即平面上 n 个点(xi,yi) i=1,…n, 寻求一个函数(曲线)y=f(x), 使 f(x) 在某种准则下与所有数据点最为接近,即曲线拟合得最好。
最常用的方法是线性最小二乘拟合 1.多项式拟合
? 对给定的数据(xj,yj ),j = 0,1,…, n;
? 选取适当阶数的多项式,如二次多项式g(x)=ax^2+bx+c;
? 使g(x)尽可能逼近(拟合)这些数据,但是不要求经过给定的数据(xj,yj ); 2.多项式拟合指令
1)多项式f(x)=a1xm+ …+amx+am+1拟合指令:
a=polyfit(x,y,m)
a :输出多项式拟合系数a[a1,a2,…,am];x,y :输出长度相同的数组;m :多项式的次数。 2)多项式在x 处的值y 的计算命令:
y=polyval(a,x )
练习4:对下面一组数据作二次多项式拟合
写出拟合命令:
plot(x,y,'k+',x,z,'r')
5
作出数据点和拟合曲线:
0.10.20.30.40.50.60.70.80.91
写出拟合的二次多项式:f (x ) =-9. 8108x +20. 1293x -0. 0317
3.可化为多项式的非线性拟和
曲线改直是工程中又一常用的判断曲线形式的方法, 许多常见的函数都可以通过适当的变换转化为线性函数。
(1)幂函数 y =ax +c
ln y -c =ln a +b ln x
b
2
c (2)指数函数 y =a b +
ln y -c =ln a =x ln b
x
+, (c (3)抛物函数 y =a x +b x
2
x ≠0 )
y -c
=ax +b x
6
练习5:完成教材P93页的习题5的第一小题。 x0=[0,300,600,1000,1500,2000]; x=0:100:2000;
y0=[0.9689,0.9322,0.8969,0.8519,0.7989,0.7491]; y=interp1(x0,y0,x,'spline'); plot(x0,y0,'k+',x,y,'r') grid on
0200400600800100012001400160018002000
7
范文四:建模 数据插值与拟合实验
数学实验与数学建模
数据插值与拟合实验
年级:20121060025 班级:电子信息科学与技术 学生姓名:吕佳琪 学号:20121060025
云南大学信息学院
一、实验目的及意义
[1] 了解插值、最小二乘拟合的基本原理
[2] 掌握用MATLAB 计算一维插值和两种二维插值的方法;
[3] 掌握用MATLAB 作最小二乘多项式拟合和曲线拟合的方法。
二、实验内容
1.针对实际问题,试建立数学模型。用MATLAB 计算一维插值和两种二维插值的方法求解;
1.用MATLAB 中的函数作一元函数的多项式拟合与曲线拟合,作出误差图;
2.用MATLAB 中的函数作二元函数的最小二乘拟合,作出误差图; 3.针对预测和确定参数的实际问题,建立数学模型,并求解。
三、实验步骤
1.开启软件平台——MATLAB ,开启MATLAB 编辑窗口; 2.根据各种数值解法步骤编写M 文件 3.保存文件并运行;
4.观察运行结果(数值或图形) ;
5.根据观察到的结果写出实验报告,并浅谈学习心得体会。
四、实验要求与任务
根据实验内容和步骤,完成以下具体实验,要求写出实验报告(实验目的→问题→数学模型→算法与编程→计算结果→分析、检验和结论→心得体会)
1.数据插值:
山区地貌:在某山区测得一些地点的高程如下表3.8。平面区域为
1200<=x>=x><><=y>=y><>
试作出该山区的地貌图和等高线图,并对几种插值方法进行比较。
表3.8 某山区高程表
y0=1200:400:3600;
z0=[1130,1250,1280,1230,1040,900,500,700; 1320,1450,1420,1400,1300,700,900,850; 1390,1500,1500,1400,900,1100,1060,950; 1500,1200,1100,1350,1450,1200,1150,1010; 1500,1200,1100,1550,1600,1550,1380,1070; 1500,1550,1600,1550,1600,1600,1600,1550; 1480,1500,1550,1510,1430,1300,1200,980];
meshc(x0,y0,z0) xlabel('x'); ylabel('y'); zlabel('z') title('原始图')
(1) 分段线性插值
[xi2,yi2]=meshgrid(1200:10:4000,1200:10:3600); zi2=interp2(x0,y0,z0,xi2,yi2,'linear'); meshc(xi2,yi2,zi2)
xlabel('x');
ylabel('y'); zlabel('z')
title('插值后')
(2) 三次样条插值
[xi,yi]=meshgrid(1200:10:4000,1200:10:3600); zi=interp2(x0,y0,z0,xi,yi,'spline'); meshc(xi,yi,zi),
title('三次样条插值表面图')
xlabel('x');ylabel('y');zlabel('z')
(3) 双三次插值
[xi,yi]=meshgrid(1200:10:4000,1200:10:3600); zi=interp2(x0,y0,z0,xi,yi,'cubic'); meshc(xi,yi,zi),
title('双三次插值后的表面图') xlabel('x轴'); ylabel('y轴'); zlabel('z轴')
(4)画等高线图:
n=9;
pcolor(x0,y0,z0) shading interp
zmax=max(max(z0));zmin=min(min(z0)); caxis([zmin,zmax]) colorbar hold on
C=contour(x0,y0,z0,n,'k:'); clabel(C)
title('未拟合图形') hold off
pcolor(xi,yi,zi) shading interp
zmax1=max(max(zi));zmin1=min(min(zi)); caxis([zmin1,zmax1]) colorbar hold on
C=contour(xi,yi,zi,n,'k:'); clabel(C) title('拟合图形') hold off
2. 曲线拟合
某年美国旧车价格的调查资料如下表所示,其中下x i 表示轿车的使用年数,
y i 表示相应的平均价格。试分析用什么形式的曲线来拟合上述的数据,并计算使用4.5年后轿车的平均价格大致为多少?
function f=E(c,x) f=c(1)*exp(-c(2)*x); 在matlab 命令窗口中输入语句:
x=[1 2 3 4 5 6 7 8 9 10];%输入x 数据
y=[2615 1943 1494 1087 765 538 484 290 226 204];%输入y 数据
c0=[1 1];%待定系数的初始值 for i=1:60%迭代系数c
c=lsqcurvefit('E' ,c0,x,y); c0=c; end
c%显示出c 的值
y1=c(1)*exp(-c(2)*x);%得出拟合函数
plot(x,y,'+',x,y1);%画出实验数据和拟合函数曲线 legend(' 实验数据' , ' 拟合曲线' )%标注图形符号意思 r=polyval(c,4.5) %求出在4.5年时的数据 Title('美国旧车数据拟合曲线')
范文五:第十讲 插值与拟合建模
第十讲 插值与拟合建模
教学内容 :插值方法与案例;拟合方法及案例
学习目标 :了解插值、拟合的基本特点;熟练掌握插值与拟合的建模案例
一、 插值问题
(一)一维 插值
1、 一维 插值的定义:已知 1+n 个节点 )
, (j
j
y x () , , 1, 0n j =)其中 j
x 互不相同,不妨设
b
x x x a n =<= 10,="" 求="" 任="" 一="" 插="" 值="">=>
*x ) (j
x ≠处的插值 *y 。
构造一个(相对简单的)函数 ) (x f y =,
通过全部节点,即:
j
j
y x f =) (, n j , , 1, 0 =
再用 ) (x f 计算插值,即 *)(*x f y =。
2、 Lagrange 插 值 :已 知 函 数
)
(x f 在
1
+n 个 点
n
x x x , , , 10 处 的 函 数 值 为
n
y y y , , , 10 . 求一 n 次多项式函数 )
(x P
n
,使其满足: n
i y x P
i i n
, , 1, 0, ) ( ==.
解决此问题的拉格朗日插值多项式公式为: ∑
=?=
n
i i
i n
y x L x P 0
) () (,其中 )
(x L
i
为 n 次多项式:
)
11101110() )(() )(() () )(() )(() (n i i i i i i i n i i i
x x x x x x x x x x x x x x x x x x x x x L
----------=
+-+- .
称为拉格朗日插值基函数。
3、 特殊情形 :(1)线性插值多项式: 10
1001
011
) (y x x x x y x x x x x L
--+
--=、
(2)抛物插值多项式:
2
120210121012002010212
)
)(() )(()
)(() )(()
)(() )(() (y x x x x x x x x y x x x x x x x x y x x x x x x x x x L ----+
----+
----=
例 1:根据下表给出的平方根值,分别用线性插值及抛物线插值计算
5
。
解 :(1)取最接近 5
=x 的两点 4
=x
, 9
1
=x
为插值节点,运用线性插值公式,
2
. 24
94539
4952) 5(51=--?
+--?
=≈L
(2) 选择与 5
=x
最接近的三点 1
=x
, 4
1
=x
, 9
2
=x
为插值节点, 根据抛物
线插值公式,有
67
. 2)
49)(19() 45)(15(3)
94)(14() 95)(15(2)
91)(41() 95)(45(1) 5(52≈----?
+----?
+----?
=≈L
4、分段性线插值 :
∑
==
n
j j j n
x l y x L
)
() (,其中:
????
?
?
???≤≤--≤≤--=+++---其它
,
0, , ) (1
1
1
1
11
j j
j j
j j
j j j
j j
x
x x
x x
x x x
x x
x
x
x x x l
5、 Matlab 作插值计算 :一维插值函数:
yi = interp1(x
, y , xi , 'method')
注意:所有的插值方法都要求 x 是单调的,并且 xi 不能够超过 x 的范围。
例 2:在 1-12的 11小时内, 每隔 1小时测量一次温度, 测得的温度依次为:5, 8,
9, 15, 25, 29, 31, 30, 22, 25, 27, 24。试估计每隔 1/10小时的温度 值。
解:hours=1:12;
temps=[5 8 9 15 25 29 31 30 22 25 27 24]; h=1:0.1:12;
t=interp1(hours,temps,h,'spline'); (直接输出数据将是很多的 )
plot(hours,temps,'+',h,t,hours,temps,'r:') % 作图 xlabel('Hour'),ylabel('Degrees Celsius’ )
(二)二维插值
1、网格节点插值法: 已知 :n m ?个节点 n
j m i z y x ij j i
, 2, 1; , , 2, 1), , , ( ==
其中 b
x x x a m =<= 21="">=>
y y y c
m =<=>=>
构造一个二函数 )
, (y x f z
=, 通过全部已知节点, 即:
ij
j i z y x f =) , ((n
j m i
, , 2, 1; , 2, 1 ==) ,再用 )
, (y x f 计算插值,即
*)
*,(*
y x f z =.
2、网格节点插值法的几种形式:
(1)最邻近点插值:
原理 :二维或高维情形的最邻近插值,与被 插值点最邻近的节点的函数值即为所求。
注意:最邻近插值一般不连续。具有连续性 的最简单的插值是分片线性插值。
(2)分片线性插值 :将四个插值点(矩 形的四个顶点)处的函数值依次简记为:
1) , (f y x f j i =, 2
1) , (f y x f j i =+
3
1
1) , (f y
x f j i =++,
4
1
) , (f y
x f j i =+
分两片的函数表达式如下:
第一片(下三角形区域) :) , (y x 满足 j
i i
i j
j y
x x x x y
y
y +---≤++) (11
插值函数为:
)
)(() )(() , (23111j i y y f f x x f f f y x f --+--+=
第二片(上三角形区域) :) , (y x 满足 j
i i
i j
j y
x x x x y
y
y +---≤++) (11
插值函数为:
)
)(() )(() , (43141i j y x f f y y f f f y x f --+--+=
注意 :) , (y x 当然应该是在插值节点所形成的矩形区域内。显然,分片线性插值函数
是连续的。
(3)双线性插值 :
特点 :双线性插值是一片一片的空间二次曲面构成。 形式 :一般形式如下:
D Cy Bx Axy y x f +++=) , (
()
)((d cy b ax ++=
)
其中有四个特定系数,利用该函数在矩形的四个顶 点(插值节点)的函数值,得到四个代数方程,正 好确定四个系数。
3、散乱节点插值法 :
已知 :1+n 个节点
)
, , (i i i
z y x
, n
i
, , 2, 1, 0 =
其中 )
, (i i
y x
互不相同。
目的 :构造一个二元函数 )
, (y x f z =, 通过全部已知 节点,即:i
i i z y x f
=) , ( (n
i
, , 1, 0 =) ,再由
)
, (y x f 计算插值,即:*)
*,(*
y x f z =.
4、 Matlab 作二维插值 :
(1) 用 MATLAB 作网格节点数据的插值:
注意:要求 x0,y0单调; x , y 可取为矩阵,或 x 取行向量, y 取为列向量, x,y 的值分别不能超出 x0,y0的范围。
例 3:测得平板表面 3*5网格点处的温度分别为:
82 81 80 82 84
79 63 61 65 81
84 84 82 85 86
试作出平板表面的温度分布曲面 z=f(x,y)的图形。
解 :1. 先在三维坐标画出原始数据,画出粗糙的温度分布曲图:
输入以下命令:
x=1:5;
y=1:3;
temps=[82 81 80 82 84;79 63 61 65 81;84 84 82 85 86];
mesh(x,y,temps)
2.以平滑数据 , 在 x 、 y 方向上每隔 0.2个单位的地方进行插值:
再输入以下命令 :
xi=1:0.2:5;
yi=1:0.2:3;
zi=interp2(x,y,temps,xi',yi,'cubic');
mesh(xi,yi,zi)
画出插值后的温度分布曲面图 .
z=interp2(x0,y0,z0,x,y,’method’)
例 4:山区地貌:在某山区测得一些地点的高程如下表。平面区域为 1200<><><><>
试作出该山区的地貌图和等高线图,并对几种插值方法进行比较。
解 :通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插值效果进行 比较。 (具体程序见课件)
(1)原始数据山区表面图 (2) 最近邻点插值山区表面图
(3) 双线性插值山区表面图 (4) 双三次插值山区表面图
(2)用 MATLAB 作散点数据的插值计算:
插值函数 griddata 格式为 :
要求 :cx 取行向量, cy 取为列向量。
二、 拟合问题
曲线拟合问题是指:已知平面上 n 个点 )
, (i i
y x , n
i
, , 2, 1 =, i
x 互不相同。寻求
函数
)
(x f ,使
)
(x f 在某种准则下与所有数据点最为接近,即曲线拟合得很好。
(一)曲线拟合最常用方法:最小二乘法 。 1、线性最小二乘法的 基本思路 : 第一步 :先选定一组函数 )
(, ), (), (21
x r x r x r m , n
m
<>
)
() () () (2211x r a x r a x r a x f m m +++= (1)
其中 m
a a a
, , , 21
为待定系数。
第二步 :确定 m
a a a
, , , 21
的准则(最小二乘准则 ) :使 n 个点 )
, (i i
y x
与曲线
)
(x f y =的距离 i
δ的平方和最小。记:
∑
∑
==-=
=
n
i i i n
i i
m y x f a a a
J 1
2
1221
]
) ([) , , , (δ
∑∑
==-=n
i i m
k i k k y x r a 1
2
1
]
) ([ (2)
问题归结为:求 m
a a a
, , , 21
使 )
, , , (21
m a a a
J 最小。
2、线性最小二乘法的求解 (1) 预备知识:
① 超定方程组:方程个数大于未知量个数的方程组。
cz =griddata(x , y , z , cx , cy , ‘method’ )
???
?
?>=+++=+++) ( 221
111212111m n y
a r a r a r y a r a r a r n m nm n n m m ,即:y a R =
其中, ????
?
?????=??????????=??
?
???????=n m
nm n n m y y y a a a r r r r r r R
11
2
1
11211
, , 。 注意 :超定方程一般是不存在解的矛盾方程组。
②如果有向量 a 使得 2
1
2211)
(i m n
i im i i y a r a r a r -+++∑
= 达到最小,则称 a 为上述超
定方程的 最小二乘解 。
(2) 线性最小二乘法的解法:曲线拟合的最小二乘法要解决的问题,实际上就是
求以下超定方程组的最小二乘解的问题:
y a R =, (3)
其中, ???
??
?????=??????????=????
??????=n m n m n m y y y a a a x r x r x r x r R
111111, , ) () () () (。
定理:当 R R T 可逆时,超定方程组(3)存在最小二乘解,且即为方程组
y
R
Ra R T
T
=
的解:y R R R a T
T
1
) (-=。
3、线性最 小二 乘拟 合 )
() () () (221
1
x r a x r a x r a x f m m +++
= 中 函数 ,
), (), ({21
x r x r
)}
(x r m 的选取:
(1)通过机理分析建立数学模型来确定 ) (x f ;
(2)将数据 ) , (i
i
y x , n i , , 2, 1 =作图,通过直观判断确定
)
(x f :
(二 ) 用 MATLAB 解拟合问题
1、用 MATLAB 作线性最小二乘拟合: (1)作多项式
1
1
21) (+-++++=m m m m
a x a x
a x
a x f 拟合 , 可利用已有程序 :
(2)
多项式在 x 处的值 y 可用以下命令计算: y = polyval(a , x )
例 1:对下面一组数据作二次多项式拟合
即要求 出二次多项式 :
322
1) (a x a x
a x f ++=中的 )
, , (321a a a A =使得 :
]) ([11
1
2
∑=-i i i
y x
f 最小。
解:解法 1.用解超定方程的方法,此时,
???
?
? ??=11 11
211
12
1x x x x R
。 1) 输入以下命令:
x=0:0.1:1;
y=[-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2]; R=[(x.^2)' x' ones(11,1)]; A=R\y'
2)计算结果: A = -9.8108 20.1293 -0.0317
0317
. 01293. 208108. 9) (2
-+-=x x x f
解法 2. 用多项式拟合的命令: 1)输入以下命令: x=0:0.1:1;
y=[-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2];
A=polyfit(x,y,2) z=polyval(A,x);
plot(x,y,'k+',x,z,'r')
% 作出数据点和拟合曲线的图形
2)计算结果:
A = -9.8108 20.1293 -0.0317
0317
. 01293. 208108. 9) (2
-+-=x x x f
三、 建模案例分析
(一) 给药方案问题
一种新药用于临床之前, 必须设计给药方案。 药物进入机体后血液输送到全身, 在这个过程中不断地被吸收、分布、代谢,最终排出体外,药物在血液中的浓度, 即单位体积血液中的药物含量,称为 血药浓度 。
一室模型:将整个机体看作一个房室,称 中心室 ,室内血药浓度是均匀的。快 速静脉注射后,浓度立即上升;然后迅速下降。当浓度太低时,达不到预期的治 疗效果;当浓度太高,又可能导致药物中毒或副作用太强。临床上,每种药物有 一个最小有效浓度 1c 和一个最大有效浓度 2c 。 设计给药方案时, 要使血药浓度 保 持在 21~c c 之间。本题设 101=c , 252=c (ug/ml)。
要设计给药方案 , 必须知道给药后血药浓度随时间变化的规律。从实验和理论 两方面着手:在实验方面 , 对某人用快速静脉注射方式一次注入该药物 300mg 后 , 在一定时刻 t(小时 ) 采集血药 , 测得血药浓度 c(ug/ml)如下表 :
解:1、 问题分析 :(1)在快速静脉注射的给药方式下,研究血药浓度(单位体积
血液中的药物含量)的变化规律。
(2) 给定药物的最小有效浓度和最大治疗浓度, 设计给药方案: 每次注射剂量多大;间隔时间多长。
实验:对血药浓度数据作拟合,符合负指数变化规律 理论:用一室模型研究血药浓度变化规律
2、模型假设:(1)机体看作一个房室,室内血药浓度均匀——一室模型; (2)药物排除速率与血药浓度成正比,比例系数 k ( > 0 ); (3)血液容积 v, t=0注射剂量 d, 血药浓度立即为 d/v。
3、模型建立 :由假设(2) ,得:
kc dt
dc -= (i ) 再由假设(3) ,得 v d c /) 0(= (ii)
解微分方程得:t
k e
v d t c ?-=)
(。
在此, d=300mg, t 及 c(t)在某些点处的值见前表,需经拟合求出参数 k 、 v 。
方法一 :用线性最小二乘拟合。
(1)先将非线性函数转化成线性函数:
)
/ln(, , ln ) /ln(ln ) (21v d a k a c y kt
v d c e
v d t c t
k =-==-=?=
?- 2
/, 121a e
d v a k a t a y =-=+=?
(2)编写 M 文件 lihe1.m 如下: d=300;
t=[0.25 0.5 1 1.5 2 3 4 6 8];
c=[19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01]; y=log(c);
a=polyfit(t,y,1) k=-a(1)
v=d/exp(a(2)) 计算结果:)
(02. 15), /1(2347. 0l v h k
==
方法二 :用非线性最小二乘拟合(略,详细过程见书本) 。
4、制定给药方案:
(1)假设:① 设每次注射剂量 D , 间隔时间 τ; ② 血药浓度 c(t) 应 21) (c t c c ≤≤;
③ 初次剂量 0D 应加大。 (2)给药方案记为:}, , {0τD D ,则:
)
( , 1220c c D c D -==νν (i )
τ
k e
c c -=21
1
2ln
1c c k
=
?τ (ii )
将 c1=10,c2=25, k=0.2347, v=15.02代入 (i)、 (ii)中,得
9. 3, 3. 225, 5. 3750===τD D
故可制定可药方案为:) (4), (225), (3750h mg D mg D ===τ;
即:首次注射 375mg ,其余每次注射 225mg ,注射的间隔时间为 4小时。
(二) 水塔流量估计问题
(一)问题的提出:
某居民区有一供居民用水的园柱形水塔,一般可以通过测量其水位来估计水 的流量,但面临的困难是,当水塔水位下降到设定的最低水位时,水泵自动启动 向水塔供水,到设定的最高水位时停止供水,这段时间无法测量水塔的水位和水 泵的供水量.通常水泵每天供水一两次,每次约两小时 .
水塔是一个高 12.2米,直径 17.4米的正园柱.按照设计,水塔水位降至约 8.2米时,水泵自动启动,水位升到约 10.8米时水泵停止工作.
表 1 是某一天的水位测量记录,试估计任何时刻(包括水泵正供水时)从水塔流 出的水流量,及一天的总用水量
表 1 水位测量记录(符号 //表示水泵启动)
(二)解题思路
1、拟合水位 ~时间函数:从测量记录看,一天有两个供水时段(以下称第 1供水时段和第 2供水时段) , 和 3个水泵不工作时段 (以下称第 1时段 t=0到 t=8.97, 第 2次时段 t=10.95到 t=20.84和第 3时段 t=23以后) .对第 1、 2时段的测量数据 直接分别作多项式拟合,得到水位函数.为使拟合曲线比较光滑,多项式次数不要太 高, 一般在 3 ~6。 由于第 3时段只有 3个测量记录, 无法对这一时段的水位作出较 好的拟合。
2、确定流量 ~时间函数:对于第 1、 2时段只需将水位函数求导数即可,对 于两个供水时段的流量,则用供水时段前后(水泵不工作时段)的流量拟合得到,并 且将拟合得到的第 2供水时段流量外推,将第 3时段流量包含在第 2供水时段内. 3、一天总用水量的估计:总用水量等于两个水泵不工作时段和两个供水时段 用水量之和,它们都可以由流量对时间的积分得到。
(三)算法设计与编程
1、拟合第 1时段的水位,并导出流量:设 t , h 为已输入的时刻和水位测量记 录(水泵启动的 4个时刻不输入) ,第 1时段各时刻的流量可如下得:
c1=polyfit(t (1:10) , h (1:10) , 3) ;
% 用 3次多项式拟合第 1时段水位, c1输出 3次多项式的系数 a1=polyder(c1) ; % a1输出多项式(系数为 c1)导数的系数
tp1=0:0.1:9;
x1=-polyval(a1, tp1) ;
% x1输出多项式 (系数为 a1) 在 tp1点的函数值 (取负后边为正值) , 即 tp1时刻的流量
于是,得到 流量函数 为:1079
.
22
7173
. 2
2356
. 0
)
(2-
+
-
=t
t
t
f
2、拟合第 2时段的水位,并导出流量:设 t , h 为已输入的时刻和水位测量记 录(水泵启动的 4个时刻不输入) ,第 2时段各时刻的流量可如下得:
c2=polyfit(t(10.9:21),h(10.9:21),3);
%用 3次多项式拟合第 2时段水位, c2输出 3次多项式的系数 a2=polyder(c2); % a2输出多项式(系数为 c2)导数的系数
tp2=10.9:0.1:21;
x2=-polyval(a2,tp2);
% x2输出多项式 (系数为 a2) 在 tp2点的函数值 (取负后边为正值) , 即 tp2时刻的流量
于是,得到 流量函数 为:
8313
. 17512. 87529. 00186. 0) (2
3--+-=t t t t f
3、拟合供水时段的流量:在第 1供水时段(t=9~11)之前(即第 1时段) 和之后 (即第 2时段) 各取几点, 其流量已经得到, 用它们拟合第 1供水时段的流量. 为 使流量函数在 t=9和 t=11连续,我们简单地只取 4个点,拟合 3次多项式(即曲线 必过这 4个点) ,实现如下:
xx1=-polyval(a1, [8 9]) ; %取第 1时段在 t=8,9的流量 xx2=-polyval(a2, [11 12]) ; %取第 2时段在 t=11,12的流量 xx12=[xx1 xx2];
c12=polyfit([8 9 11 12], xx12, 3) ; %拟合 3次多项式 tp12=9:0.1:11;
x12=polyval(c12, tp12) ; % x12输出第 1供水时段各时刻的流量
在第 2供水时段之前取 t=20, 20.8两点的流水量,在该时刻之后(第 3时段) 仅有 3个水位记录,我们用差分得到流量,然后用这 4个数值拟合第 2供水时段的流 量如下:
dt3=diff(t(22:24) ) ; %最后 3个时刻的两两之差 dh3=diff(h(22:24) ) ; %最后 3个水位的两两之差 dht3=-dh3./dt3; %t(22)和 t(23)的流量 t3=[20 20.8 t(22) t(23)];
xx3=[-polyval(a2, t3(1:2) , dht3)]; %取 t3各时刻的流量 c3=polyfit(t3, xx3, 3) ; %拟合 3次多项式 t3=20.8:0.1:24; x3=polyval(c3, tp3) ; % x3输出第 2供水时段(外推至 t=24)各时刻的流量 拟合的流量函数为:
8283. 913077. 71405. 0) (2
-+-=t t t f 。
4、 一天总用水量的估计:第 1、 2时段和第 1、 2供水时段流量的积分之和,就 是一天总用水量.虽然诸时段的流量已表为多项式函数,积分可以解析地算出,这里 仍用数值积分计算如下:
y1=0.1*trapz(x1); %第 1时段用水量(仍按高度计) , 0.1为积分步长
y2=0.1*trapz(x2); %第 2时段用水量
y12=0.1*trapz(x12); %第 1供水时段用水量 y3=0.1*trapz(x3); %第 2供水时段用水量
y=(y1+y2+y12+y3) *237.8*0.01; %一天总用水量 L m 3310=
计算结果:y1=146.2, y2=266.8, y12=47.4, y3=77.3, y=1250.4
5、 流量及总用水量的检验:计算出的各时刻的流量可用水位记录的数值微分来 检验.用水量 y1可用第 1时段水位测量记录中下降高度 968-822=146来检验,类似 地, y2用 1082-822=260检验.
供水时段流量的一种检验方法如下:供水时段的用水量加上水位上升值 260是 该时段泵入的水量,除以时段长度得到水泵的功率(单位时间泵入的水量) ,而两个 供水时段水泵的功率应大致相等.第 1、 2时段水泵的功率可计算如下:
p1=(y12+260)/2; %第 1供水时段水泵的功率(水量仍以高度计)
tp4=20.8:0.1:23; xp2=polyval(c3, tp4) ; % xp2输出第 2供水时段各时刻的流量
p2=(0.1*trapz(xp2)+260)/2.2; %第 2供水时段水泵的功率(水量仍以高度计)
计算结果 :p1=154.5 , p2=140.1
(四)计算结果
流量函数为:
??????
?≤≤-+-
<><>
<>
2421 8283. 91 3077. 7 1405.
021
11 8313. 1 7512. 8 7529. 0 0186. 0 119 078. 355 5879. 73 7207. 3 90 1079. 22 7173. 2 2356. 0) (22
3
2
2
t t t t t t t t t t t t t t f
流量曲线见图 :
n = (3,4) n = (5,6)