教案课件数值计算一维

1、钮海钮海考试提要考试提要晚上行政搂后面的工程设计中心工程设计中心，光纤网带上学生证，办卡后1元1小时。软件网址：202.115.42.124/home/model用户名：model联系Email：Tel:89859857上机上机组队组队n编队配组：编队配组：建模，计算，写作建模，计算，写作n自己组队：自己组队：3人人/队队n自愿参加：自愿参加：300元元/队队 报名费报名费请发请发email：写明：团队每个成员的写明：团队每个成员的 姓名，学院，性别，手机姓名，学院，性别，手机Excel与与Matlab数据传递数据传递位于位于matlabtoolboxexlink目录下目录下前言前言n内容：内

2、容： 一维问题一维问题 高维问题高维问题n插值插值n拟合拟合n插值插值n回归回归前言前言n层次：层次： 插值与拟合插值与拟合理论理论 Matlab实现实现 应用应用实例实例前言前言问题问题n测试数据测试数据自变量自变量数据数据因变量因变量数据数据n理论研究理论研究离散函数离散函数因变量因变量f(自变量自变量)连续函数连续函数分类分类n一个自变量一个自变量一维问题一维问题n多个自变量多个自变量多维问题多维问题方法方法n1.测试数据为测试数据为精确精确数据或者数据或者误差可以忽略误差可以忽略 方法：方法：插值插值法法n2.测试数据含测试数据含误差误差，必须考虑，必须考虑误差误差 方法：方法：拟合拟

3、合法或法或回归回归一维问题一维问题插值插值一维问题一维问题-插值插值xyy1y2y3y4x1x2x3x4xy？一维问题一维问题-插值插值n插值：插值：xX1 X2 Xn yy1 y2 yn 问题问题n要求：对任意自变量要求：对任意自变量x1xxn，求，求x对应的对应的y值值插值点插值点插值插值n重点：重点：找函数关系找函数关系f，使得，使得yf(x)转换转换一维问题一维问题-插值插值n首要问题：首要问题： f的选取的选取？多项式函数多项式函数1110( )mmmma xaxafaxx满足：满足：yif(xi)，i1n小要求：小要求：多项式函数多项式函数唯一唯一！数学模型数学模型一维问题一维问题

4、-插值插值n求解：求解：11101.,.mmmmaaayxixiinaxiiyif(xi)，i1n线性方程组线性方程组未知变量个数：未知变量个数：m+1方程个数：方程个数：n唯一唯一m+1=nm=n-1一维问题一维问题-插值插值n解：解：101111. 112. 1.1nnnnxyxynxnaa一维问题一维问题-插值插值n模型分析：模型分析：Runge现象现象 测试点越多，多项式次数越高测试点越多，多项式次数越高失真失真21,(1) 5,5,11xnyx 10( )Px-505-1.5-1-0.500.511.52y=1/(1+x2)n=2n=4n=6n=8n=10一维问题一维问题-插值插值n

5、失真失真解决办法：解决办法：每次采取数据点比较少每次采取数据点比较少 如如2个；个；4个个n问题：原先大量数据？问题：原先大量数据？n方法：分组（每组方法：分组（每组2个；个；4个个递进进行递进进行）分段线性插值分段线性插值linear分段三次插值分段三次插值cubic一维问题一维问题-插值插值n再处理！（光滑）再处理！（光滑）n可导可导 要求要求：（：（见后注见后注） 1.插值函数插值函数f（分段函数分段函数）在整个在整个x1,xn具有二阶具有二阶连续导数；连续导数； 2.在每个小区间在每个小区间xi-1,xi是三次多项式是三次多项式n新问题不够完美，新问题不够完美，光滑光滑？一维问题一维问

6、题-插值插值n数学模型：数学模型：yif（xi），），i1nf（x），），f （x），），f （x）在）在xi连续，连续，i2n-1进一步表示（省略）进一步表示（省略）n求解求解：省略省略一维问题一维问题-插值插值n后注：后注：取自工程中利用取自工程中利用弹性竹条弹性竹条（样条）（样条）进行插值的进行插值的方法方法n因此因此该插值称为该插值称为样条样条spline插值方法插值方法实现实现n编程？编程？nMatlab实现实现：插值函数插值函数（程序）（程序）插值插值=interp1(自变量数据自变量数据,因变量试数据因变量试数据,插值点插值点,方法方法)方法方法:1.省却省却linear分段线分

7、段线性性 2.cubic分段三次分段三次 3.spline样条样条实验实验问题问题21,(1) 5,5,10 xnyx x-5 5y 0.038 0.038分别用分别用分段线性、分段三次、样条分段线性、分段三次、样条插值插值与原函数比较与原函数比较实验实验脚本文件脚本文件n%各类一维分段插值比较脚本文件脚本文件（文件名（文件名cz.mcz.m）nx=linspace(-5,5,9);%原始自变量ny=1./(1+x.2);%原始因变量nx_o=linspace(-5,5,50);ny_o=1./(1+x_o.2);%原始函数nx_in_l=linspace(-5,5,50);ny_in_l=i

8、nterp1(x,y,x_in_l);%线性插值函数nsubplot(3,1,1)%子图nplot(x_o,y_o,x_in_l,y_in_l)ntitle(linear)nx_in_c=x_in_l;ny_in_c=interp1(x,y,x_in_c,cubic);%三次插值函数nsubplot(3,1,2)nplot(x_o,y_o,x_in_c,y_in_c)ntitle(cubic)nx_in_s=x_in_l;ny_in_s=interp1(x,y,x_in_s,spline);%样条插值函数nsubplot(3,1,3)nplot(x_o,y_o,x_in_s,y_in_s)nt

9、itle(spline)实验实验运行运行n在在matlab工作区输入工作区输入cz-5-4-3-2-101234500.51linear-5-4-3-2-101234500.51cubic-5-4-3-2-101234500.51spline保形保形不光滑不光滑光滑光滑不保形不保形插值方法选择插值方法选择n插值函数插值函数插值插值=interp1(自变量数据自变量数据,因变量试数据因变量试数据,插值点插值点,方法方法)方法方法:1.省却省却linear分段线分段线性性 2.cubic分段三次分段三次 3.spline样条样条要求保形要求保形快速粗估快速粗估要求光滑要求光滑建模实例建模实例n估计

10、水塔的水流量估计水塔的水流量(AMCM92A)n美国某州的各用水管理机构要求各社区提供以每小时多少加仑计的用水美国某州的各用水管理机构要求各社区提供以每小时多少加仑计的用水率以及每天所用的总水量。但许多社区并没有测量流入或流出当地水塔率以及每天所用的总水量。但许多社区并没有测量流入或流出当地水塔的水量的设备，他们只能代之以每小时测量水塔中的水位，其精度在的水量的设备，他们只能代之以每小时测量水塔中的水位，其精度在5%以内。更为重要的是，无论什么时候，只要水塔中的水位下降到某一最以内。更为重要的是，无论什么时候，只要水塔中的水位下降到某一最低水位低水位L时，水泵就启动向水塔重新充水直至某一最高水

11、位时，水泵就启动向水塔重新充水直至某一最高水位H，但也无法，但也无法得到水泵的供水量的测量数据。因此，在水泵正在工作时，人们不容易得到水泵的供水量的测量数据。因此，在水泵正在工作时，人们不容易建立水塔中的水位与水泵工作时的用水量之间的关系。水泵每天向水塔建立水塔中的水位与水泵工作时的用水量之间的关系。水泵每天向水塔充水两次，每次约二小时。充水两次，每次约二小时。n试估计在任何时刻，甚至包括水泵正在工作的时间内，水从水塔流出的试估计在任何时刻，甚至包括水泵正在工作的时间内，水从水塔流出的流量流量f (t)，并估计一天的总用水量和水泵的工作功率。表，并估计一天的总用水量和水泵的工作功率。表1给出了

12、某个真给出了某个真实小镇某一天的真实数据。实小镇某一天的真实数据。建模实例建模实例时间时间（秒）（秒）水位水位(0.01英尺英尺)时间时间（秒）（秒）水位水位(0.01英尺英尺)时间时间（秒）（秒）水位水位(0.01英尺英尺)0317535932水泵工作6853528423316311039332水泵工作7185427676635305439435355075021269710619299443318344579254水泵工作13937294746636335082649水泵工作1792128924995332608596834752124028505393631678995333972522

13、32797572543087932703340285432752605743012322842697645542927表表1 某小镇某天的水塔水位某小镇某天的水塔水位 建模实例建模实例n表表1给出了从第一次测量开始的以秒为单位的给出了从第一次测量开始的以秒为单位的时刻，以及该时刻的高度单位为百分之一英尺时刻，以及该时刻的高度单位为百分之一英尺的水塔中水位的测量值，例如，的水塔中水位的测量值，例如，3316秒后，水秒后，水塔中的水位达到塔中的水位达到31.10英尺。水塔是一个垂直英尺。水塔是一个垂直圆形柱体，高为圆形柱体，高为40英尺，直径为英尺，直径为57英尺。通常英尺。通常当水塔的水位降至约

14、当水塔的水位降至约27.00英尺时水泵开始向英尺时水泵开始向水塔充水，而当水塔的水位生至约水塔充水，而当水塔的水位生至约35.50英尺英尺时水泵停止工作。时水泵停止工作。建模实例建模实例nt0=0,3316,6635,10619,13937,17921,21240,25223,28543,32284,39435,43318,46636,49953,. 53936,57254,60574,64554,68535,71854,75021,85968,89953,93270;nL0=3175,3110,3054,2994,2947,2892,2850,2797,2752,2697,3550,3445

15、,3350,3260,3167,3087,. 3012,2927,2842,2767,2697,3475,3397,3340;nt=linspace(0,86400,200);nL=interp1(t0,L0,t,clubic);nplot(t0,L0,*,t,L)建模实例建模实例012345678910 x 104260028003000320034003600spline012345678910 x 104260028003000320034003600cubic一维问题一维问题拟合拟合一维问题一维问题-拟合拟合xyy1y2y3y4x1x2x3x4？323210a xa xa xa3223

16、yxxx3221yxxx一维问题一维问题-拟合拟合n拟合：拟合：xX1 X2 Xn yy1 y2 yn 问题问题参数参数n要求：要求：理论推导理论推导或或经验估计经验估计函数关系为函数关系为yf(x,t)n重点：重点：找最合适的找最合适的t，使理论和实际值，使理论和实际值之间的之间的总体误差最小总体误差最小理论理论实际实际多项式系数多项式系数一维问题一维问题-拟合拟合xyy1y2y3y4x1x2x3x4一维问题一维问题-拟合拟合n首要问题：首要问题： f的选取的选取？根据问题性质根据问题性质( , )f x t误差：误差：1.单点误差：单点误差：Ciyif(xi,t)，i1n2.总体误差：总体

17、误差： 或或 等等要求：要求：总体误差总体误差最小最小！数学模型数学模型1|niiC21niiC21min( , )niiiyf x tt最小二乘法最小二乘法n求解：求解：二次无约束规划问题（二次无约束规划问题（多元函数极值多元函数极值）一维问题一维问题-拟合拟合21min( , )niiiyf x tt理论上：理论上：找驻点找驻点数学模型数学模型实现实现多项式拟合多项式拟合n编程？编程？nMatlab实现实现：1.多项式多项式拟合拟合（程序）（程序）拟合多项式系数向量（高拟合多项式系数向量（高低）低）=polyfit(自变量数据自变量数据,因变量试数据因变量试数据,拟合多项式次数拟合多项式次

18、数)拟合多项式次数拟合多项式次数1，线性拟合（回归），线性拟合（回归）实验实验问题问题t00.30.81.11.62.3 y0.50.821.141.251.351.40 用用多项式多项式拟合拟合实验实验上机观察上机观察nt = 0，.3，.8 ，1.1，1.6 ，2.3; ny = 0.5，0.82，1.14，1.25，1.35，1.40; nplot(t,y,o) ngrid on二次二次实验实验分析分析n多项式回归多项式回归n由图可以看出应该可以用二次多项式来由图可以看出应该可以用二次多项式来表达表达：y=a2*t2+a1*t +a0实验实验操作操作napolyfit(t,y,2)na

19、-0.2387 0.9191 0.5318 ntf=linspace(0,2.3,50);nyf=polyval(a,x);nplot(t,y,o,tf,yf)ngrid on实验实验操作操作n结果令人失望，结果令人失望，但我们可以增加但我们可以增加阶数来提高精确阶数来提高精确度，但更明智的度，但更明智的选择是用别的方选择是用别的方法法 实验实验改进改进n用指数关系表达用指数关系表达： y=a0+a1*exp(-t)+a2*exp(-t)2 实验实验操作操作nxexp(-t);na polyfit(x,y,2)na 0.4097 0.8988 1.3974ntf=linspace(0,2.3,

20、50);nyf=polyval(a,x);nplot(t,y,o,tf,yf)ngrid on实验实验操作操作n看起来是不是看起来是不是好多了好多了 实验实验1790 1800 1810 1820 1830 1950 1960 1970 1980 3.9 5.3 7.2 9.6 12.9 150.7 179.3 204.0 226.5例：美国人口数据（单位例：美国人口数据（单位1.0e+06 ）cdate pop 用多项式拟合用多项式拟合实验实验np = polyfit(cdate,pop,4) nWarning: Matrix is close to singular or badly sc

21、aled. Results may be inaccurate. RCOND = 5.429790e20 p = 1.0e+05 * 0.0000 0.0000 0.0000 0.0126 6.0020 实验实验n产生警告的原因是计算中的产生警告的原因是计算中的cdata值太大，在值太大，在计算中的计算中的Vandermonde行列式使变换产生了行列式使变换产生了问题，解决的方法之一是使数据标准化问题，解决的方法之一是使数据标准化n标准化方法：标准化方法：data-Mean (data)/Std(data)实验实验n数据的标准化数据的标准化是对数据进行缩放，以使以后的是对数据进行缩放，以使以后

22、的计算能更加精确，一种方法是使之成为计算能更加精确，一种方法是使之成为0均值：均值：nmatlab实现方法：实现方法：nsdate = (cdate mean(cdate)./std(cdate) 实验实验np = polyfit(sdate,pop,4) p = 0.7047 0.9210 23.4706 73.8598 62.2285npop4 = polyval(p,sdate); nplot(cdate,pop4,cdate,pop,+), grid on 实验实验 1.数据类型的一致化数据类型的一致化 极大型指标极大型指标: :总是期望指标的取值越大越好；总是期望指标的取值越大越好；

23、 极小型指标极小型指标: :总是期望指标的取值越小越好；总是期望指标的取值越小越好； 中间型指标中间型指标: :总是期望指标的取值既不要太大，也不要总是期望指标的取值既不要太大，也不要太小为好，即取适当的中间值为最好太小为好，即取适当的中间值为最好; ; 区间型指标区间型指标: :总是期望指标的取值最好是落在某一个确总是期望指标的取值最好是落在某一个确定的区间内为最好。定的区间内为最好。 数据的规范化处理数据的规范化处理 （1） 极小型指标） 极小型指标: 对于某个极小型指标对于某个极小型指标x， 则通过变换， 则通过变换1(0)xxx , ,或变换或变换xMx , ,其中其中M为指标为指标x

24、的可能的可能取值的最大值，即可将指标取值的最大值，即可将指标x极大化。极大化。 数据数据类型的一致化类型的一致化 （2）中中间间型型指指标标: 对对于于某某个个中中间间型型指指标标x，则则通通过过变变换换 2()1,()22()1,()2xmmxMmMmxMxMmxMMm 其其中中M和和m分分别别为为指指标标x的的可可能能取取值值的的最最大大值值和和最最小小值值，即即可可将将中中间间型型指指标标x极极大大化化。 指标类型的一致化指标类型的一致化 （3）区区间间型型指指标标 对对于于某某个个区区间间型型指指标标x，则则通通过过变变换换 1,1,1,axxacxaxbx bxbc 其其中中 , a

25、 b为为指指标标x的的最最佳佳稳稳定定的的区区间间，max,ca m Mb，M和和m分分别别为为指指标标x的的可可能能取取值值的的最最大大值值和和最最小小值值。即即可可将将区区间间型型指指标标x极极大大化化。 2. 评价指标的无量纲化评价指标的无量纲化 指标的规范化处理指标的规范化处理 在在实实际际中中的的评评价价指指标标12,(1)mx xxm之之间间，往往往往都都存存在在着着各各自自不不同同的的单单位位和和数数量量级级，使使得得这这些些指指标标之之间间存存在在着着不不可可公公度度性性，这这就就为为综综合合评评价价带带来来了了困困难难，尤尤其其是是为为综综合合评评价价指指标标建建立立和和依依

26、据据这这个个指指标标的的大大小小排排序序产产生生不不合合理理性性。 如果不对这些指标作相应的无量纲处理，则在综合评如果不对这些指标作相应的无量纲处理，则在综合评价过程中就会出价过程中就会出“大数吃小数大数吃小数”的错误结果，从而导致最的错误结果，从而导致最后得到错误的评价结论。后得到错误的评价结论。 无量纲化处理又称为指标数据的无量纲化处理又称为指标数据的标准化标准化,或或规范化规范化处理。处理。 常用方法常用方法:标准差方法、极值差方法和功效系数方法等。标准差方法、极值差方法和功效系数方法等。 指标的无量纲化指标的无量纲化 假假设设m个个评评价价指指标标12,mx xx，在在此此不不妨妨假假

27、设设已已进进行行了了类类 型型 的的 一一 致致 化化 处处 理理 ， 并并 都都 有有n组组 样样 本本 观观 测测 值值(1,2, ;1,2, )ijx in jm，则则将将其其作作无无量量纲纲化化处处理理。 （1）标标准准差差方方法法: 令令ijjijjxxxs (1,2, ;1,2, )in jm， 其其中中1221111,() (1,2, )nnjijjijjiixx sxxjmnn。 显显然然指指标标(1,2, ;1,2, )ijx in jm的的均均值值和和均均方方差差分分别别为为0 0和和1 1，即即0,1ijx 是是无无量量纲纲的的指指标标，称称之之为为ijx的的标标准准观观

28、测测值值。 评价指标的无量纲化评价指标的无量纲化 （2）极极值值差差方方法法: 令令ijjijjjxmxMm (1,2, ;1,2, )in jm， 其其中中11max ,min (1,2, )jijjiji ni nMxmxjm 。则则0,1ijx 是是无无量量纲纲的的指指标标观观测测值值。 （3）功功效效系系数数法法: 令令ijjijjjxmxcdMm (1,2, , ;1,2, , )in jm， 其其中中, c d均均为为确确定定的的常常数数。c表表示示“平平移移量量” ，d表表示示“旋旋转转量量” ，即即表表示示“放放大大”或或“缩缩小小”倍倍数数，则则 ,ijxc c d 。 譬譬

29、如如若若取取60,40cd，则则60,100ijx 。 nmatlab中的归一化处理有三种方法中的归一化处理有三种方法n1. premnmx、postmnmx、tramnmxn2. restd、poststd、trastd实现实现曲线拟合曲线拟合nMatlab实现实现：1.最小二乘最小二乘拟合拟合程序程序参数参数=lsqnonlin(误差函数误差函数, 参数估计值参数估计值)误差函数文件误差函数文件function c=myfun(参数参数)自变量数据自变量数据=x1,xn;因变量数据因变量数据=y1,yn;c=f(自变量数据自变量数据,参数参数)-因变量数据因变量数据;拟合函数数学表达式（拟

30、合函数数学表达式（向量形式向量形式）建模实例建模实例 t (h) 0.25 0.5 1 1.5 2 3 4 6 8c ( g/ml) 19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01已知一室模型快速静脉注射下的血药浓度数据已知一室模型快速静脉注射下的血药浓度数据(t=0注射注射300mg)求血药浓度随时间的变化规律求血药浓度随时间的变化规律c(t).建模实例建模实例n解：解：0( )ktc tc e由机理由机理微分方程微分方程其中，其中，c0，k为待定参数为待定参数建模实例建模实例n求解：求解：n1.建立误差函数文件：文件名为建立误差函数文件：

31、文件名为nd_f.mfunction wc=nd_f(x)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;wc= x(1)*exp(-x(2)*t)c;0ktc e其中：其中： c0 x(1), k x(2);拟合函数为拟合函数为建模实例建模实例n2.在工作区输入：在工作区输入：nx=lsqnonlin(nd_f, 1,1)nx= 20.2413 0.2420即即0.24( )20.24tc te建模实例建模实例n3.在工作区输入：在工作区输入：nt=0.25 0.5 1 1.5 2

32、3 4 6 8;nc=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01;ntf=linspace(0,10,50);ncf=x(1)*exp(-x(2)*tf);nplot(t,c,o,tf,cf)建模实例建模实例0123456789100510152025实现实现曲线拟合曲线拟合nMatlab实现实现：2.曲线曲线拟合拟合程序程序参数参数= lsqcurvefit(拟合函数拟合函数, 参数估计值参数估计值,自变量数据自变量数据,因变量数据因变量数据)拟合函数文件拟合函数文件function y=myfun(参数参数,自变量自变量)y=f(自变

33、量自变量,参数参数);拟合函数数学表达式（拟合函数数学表达式（向量形式向量形式）nlsqcurvefit( (参数参数,自变量自变量) 拟合函数拟合函数(自变自变量量,参数参数), 参数估计值参数估计值,自变量数据自变量数据,因变量数据因变量数据)建模实例建模实例n求解：求解：n1.建立误差函数文件：文件名为建立误差函数文件：文件名为nd1_f.mfunction c=nd1_f(x,t)c= x(1)*exp(-x(2)*t);0ktc e其中：其中： c0 x(1), k x(2);拟合函数为拟合函数为建模实例建模实例n在工作区输入：在工作区输入：nt=0.25 0.5 1 1.5 2 3

34、 4 6 8;nc=19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01;nx= lsqcurvefit(nd1_f,1,1,t,c)nx= 20.2413 0.24200ktc e拟合好坏的依据拟合好坏的依据n在由经验估计的拟合函数中在由经验估计的拟合函数中n估计拟合函数好坏的评判标准有：估计拟合函数好坏的评判标准有： 1.比较总体误差大小比较总体误差大小n总体误差越小越好总体误差越小越好 2.观察误差观察误差余量分析余量分析n理想的误差应该是理想的误差应该是随机随机的的实验实验1790 1800 1810 1820 1830 1950 1960

35、 1970 1980 3.9 5.3 7.2 9.6 12.9 150.7 179.3 204.0 226.5例：美国人口数据（单位例：美国人口数据（单位1.0e+06 ）cdate pop 用多项式拟合用多项式拟合实验实验n线性拟合线性拟合np1 = polyfit(sdate,pop,1);npop1 = polyval(p1,sdate);nplot(cdate,pop1,cdate,pop,+) 实验实验实验实验n余量分析余量分析 nres1 = pop pop1;nFiguren plot(cdate,res1,+) 实验实验n抛物线拟合抛物线拟合np2 = polyfit(sdat

36、e,pop,2);npop2= polyval(p2,sdate);nplot(cdate,pop2,cdate,pop,+) 实验实验实验实验n余量分析余量分析 nres2 = pop pop2;nFiguren plot(cdate,res2,+) 实验实验n四次拟合四次拟合np4 = polyfit(sdate,pop,4);npop4= polyval(p4,sdate);nplot(cdate,pop4,cdate,pop,+) 实验实验实验实验n余量分析余量分析 nres4 = pop pop4;nFiguren plot(cdate,res4,+) 可以看出，多项式拟合即使可以看出，多