实验六、拟合与插值问题

1、2021/6/161数学模型数学模型数学实验数学实验 六六一、实验项目名称一、实验项目名称 传染病模型传染病模型 插值与拟合插值与拟合2021/6/162进一步巩固、加强进一步巩固、加强Matlab的应用能力、的应用能力、学会用学会用MATLAB软件进行数据拟合软件进行数据拟合了解在最小二乘意义下数据拟合的理了解在最小二乘意义下数据拟合的理论和方法论和方法.通过对实际问题的分析和研究，初步通过对实际问题的分析和研究，初步掌握建立数据拟合数学模型的方法掌握建立数据拟合数学模型的方法插值的基本原理插值的基本原理二、实验目的二、实验目的2021/6/163三、实验内容与步三、实验内容与步骤骤1 1、

2、建模实例：传染、建模实例：传染病模型病模型2 2、MATLABMATLAB求函求函数的拟合与插值数的拟合与插值2021/6/164 目目 的的1、传染病模型、传染病模型 描述传染病的传播过程描述传染病的传播过程 分析受感染人数的变化规律分析受感染人数的变化规律 预报传染病高潮到来的时刻预报传染病高潮到来的时刻 预防传染病蔓延的手段预防传染病蔓延的手段 按照传播过程的一般规律，按照传播过程的一般规律，用机理分析方法建立模型用机理分析方法建立模型2021/6/1651.1问题重述问题重述n问题: 有一种传染病（如SARS、甲型H1N1）正在流行。现在希望建立适当的数学模型，利用已经掌握的一些数据资

3、料对该传染病进行有效地研究，以期对其传播蔓延进行必要的控制，减少人民生命财产的损失。考虑如下的几个问题，建立适当的数学模型，并进行一定的比较分析和评价展望。n1、不考虑环境的限制，设单位时间内感染人数的增长率是常数，建立模型求t时刻的感染人数。n2、假设环境条件下所允许的最大可感染人数为 。单位时间内感染人数的增长率是感染人数的线性函数，最大感染时的增长率为零。建立模型求t时刻的感染人数。mx2021/6/1663、现有卫生防疫部门采集到的某地区一定时间内一定间隔区间的感染人数数据（见下表），利用该数据确定上述两个模型中的相关参数，并将它们的预测值与实际数据进行比较分析（计算仿真偏差）并对两个

4、模型进行适当的评价。（注：该问题中，设最大可感染人数为2000人）2021/6/167实验问题实验问题2021/6/168模型（一）模型（一）模型（二）模型（二）Matlab-Matlab-微分方程的求解微分方程的求解已解决已解决今天我们研究问今天我们研究问题题4 4的方法的方法拟合拟合插值插值2021/6/169据人口统计年鉴，知我国从据人口统计年鉴，知我国从19491949年至年至19941994年人年人口数据资料如下：口数据资料如下： ( (人口数单位为：百万人口数单位为：百万) )（1 1）在直角坐标系上作出人口数的图象。）在直角坐标系上作出人口数的图象。（2 2）建立人口数与年份的函

5、数关系，并估算）建立人口数与年份的函数关系，并估算19991999年年的人口数。的人口数。拟拟 合合 问问 题题 实实 例例 1年份年份19491954 1959 1964 1969人口数人口数 541.67602.66 672.09 704.99 806.71 年份年份 1974 1979 1984 1989 1994人口数人口数 908.59 975.42 1034.751106.761176.74 2021/6/1610如何确定如何确定21axay线性模型线性模型21,aa2021/6/16111 曲线拟合问题的提法曲线拟合问题的提法: 已知一组（二维）数据，即平面上的已知一组（二维）数

6、据，即平面上的n 个点个点),(iiyx， ixni, 2 , 1L互不相同，寻求一个函数（曲线）互不相同，寻求一个函数（曲线）)(xfy ，使使)(xf在某种准则下与所有数据点最为接近，即曲线拟合在某种准则下与所有数据点最为接近，即曲线拟合得最好，如图得最好，如图: xy0+i ),(iiyx)(xfy 一、曲线拟合一、曲线拟合 niiiniiyxf1212)( 确定确定f(x)使得使得 达到最小达到最小 最小二乘准则最小二乘准则 2021/6/1612. 用什么样的曲线拟合已知数据用什么样的曲线拟合已知数据?常用的曲线函数系类型：常用的曲线函数系类型：画图观察；画图观察；理论分析理论分析x

7、aeay21 指数曲线：指数曲线： 21axay 双曲线（一支双曲线（一支):): 11 mmmaxaxayL多项式：多项式： 21axay直线：直线： 2021/6/1613线性最小二乘拟合线性最小二乘拟合 f(x)=a1r1(x)+ +amrm(x)中中函数函数 r1 1( (x),), ,rm( (x)的选取的选取 1. 1. 通过机理分析建立数学模型来确定通过机理分析建立数学模型来确定 f( (x) )；+f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx 2. 2. 将数据将数据 (xi,yi) i=1, ,n 作图，通过

8、直观判断确定作图，通过直观判断确定 f(x)：2021/6/1614 拟合函数组中系数的确定拟合函数组中系数的确定达达到到最最小小。 niiiniiyaxaaaJ12211221),( 21,aa,)(21为例为例以以axaxf 使使得得，即即确确定定21aa 0) 20) 2121121niiiiniiiyaxaxyaxa2021/6/1615用用MATLAB作线性最小二乘拟合作线性最小二乘拟合1. 1. 作多项式作多项式f(x)=a1xm+ +amx+am+1拟合拟合, ,可利用已有程序可利用已有程序:a=polyfit(x,y,m)2.2.多项式在多项式在x处的值处的值y可用以下命令计算

9、：可用以下命令计算： y=polyval（a，x）输出拟合多项式系数输出拟合多项式系数a=a1, ,am , am+1 (数组）数组）输入同长度输入同长度的数组的数组x，y拟合多项拟合多项式次数式次数2021/6/1616人口预测线性模型人口预测线性模型对于开始提出的实验问题对于开始提出的实验问题, 代如数据，计算得代如数据，计算得27754,1521 aa从而得到人口数与年份的函数关系为从而得到人口数与年份的函数关系为2775415 xy把把x=1999代如，估算出代如，估算出1999年的人口数为年的人口数为 y=1252.1（百万）（百万）12.52亿亿1999年实际人口数量为年实际人口数

10、量为.亿。亿。线性预测模型线性预测模型2021/6/1617syms xx=1949:5:1994;y=541.67 602.66 672.09 704.99 806.71 908.59 975.42 1034.75 1106.76 1176.74;plot(x,y,r*)z=polyfit(x,y,1)y1=polyval(z,x)hold onplot(x,y1,b+-)19451950195519601965197019751980198519901995500600700800900100011001200人口模型的解人口模型的解2021/6/1618nsyms xnx=0:1:14n

11、y=39 53 72 96 129 171 232 314 386 502 629 760 920 1065 1232;nplot(x,y,r*)nz=polyfit(x,y,2)ny1=polyval(z,x)nhold onnplot(x,y1,b+-)传染病模型传染病模型3 3的解的解0246810121402004006008001000120014002021/6/1619二、插二、插 值值 2. 面面 对对 一一 个个 实实 际际 问问 题，应题，应 该该 用用 插插 值，还值，还 是是 拟拟 合。合。1. 插值的基本原理；三种插值方法：拉格朗日插值的基本原理；三种插值方法：拉格朗

12、日插值，分段线性插值，三次样条插值。插值，分段线性插值，三次样条插值。2021/6/1620插插 值值 问问 题题 实实 例例 1 1xy机翼下轮廓线2021/6/1621插插 值值 问问 题题 的的 提提 法法已知已知 n+1个节点个节点, 1 , 0(),(njyxjjL其中其中jx互不相同，不妨设互不相同，不妨设),10bxxxanL求任一插值点求任一插值点)(*jxx 处的插值处的插值.*y0 x1xnx0y1y节点可视为由节点可视为由)(xgy 产生产生,g表达式复杂表达式复杂或无封闭形式或无封闭形式,，或未知或未知.。*x*y2021/6/16220 x1xnx0y1y求求 解解

13、插插 值值 问问 题题 的的 基基 本本 思思 路路 构造一个构造一个(相对简单的相对简单的)函数函数),(xfy 通过全部节点通过全部节点, 即即), 1 ,0()(njyxfjjL再用再用)(xf计算插值，即计算插值，即).(*xfy *x*y2021/6/1623插值的基本原理插值的基本原理选讲选讲PPT（24-30）2021/6/16241.1.拉格朗日(Lagrange)多项式插值1.1 1.1 插值多项式插值多项式) 1 ()(0111axaxaxaxLnnnnnLnnnnnnnnyyYaaAxxxxXLLL001100,11)2(,0)det(X有唯一解有唯一解), 1 , 0(

14、)(njyxLjjnL)2(YXA 2021/6/16251.1.拉格朗日(Lagrange)多项式插值1.2 拉格朗插值多项式nixxxxxxxxxxxxxxxxxlniiiiiiniiiLLLLL1 , 0,)()()()()()()(110110)3()()(0 xlyxLiniinjjnjiyxLjijixl)(,0, 1)(又（2）有唯一解，故（3）与（1）相同。2021/6/16261.1.拉格朗日(Lagrange)多项式插值1.3 1.3 误差估计误差估计),(),()!1()()()()(0)1(baxxngxLxgxRnjjnnn1)1()(nnMgnjjnnxxnMxR0

15、1)!1()()(xRnn)(xRgn光滑)(xRxxnj接近2021/6/16271.1.拉格朗日(Lagrange)多项式插值1.4 1.4 例例 将0,/2 n等分，用g(x)=cos(x)产生n+1个节点，作Ln(x)（取n=1,2) ,计算cos(/6), 估计误差。解解: n=1, (x0,y0)=(0,1), (x1,y1)=(/2,0), L1(x)=y0l0+y1l1=1-2x/, cos(cos( /6)=0.6667/6)=0.6667n=2, (x0,y0)=(0,1), (x1,y1)=(/4,0.7071),(x2,y2)=(/2,0),L2(x)=y0l0+y1l

16、1+y2l2=8(x-/4)(x-/2)/2-16x(x-/2)0.7071/2cos(cos( /6)=L/6)=L2 2( ( /6)=0.8508 /6)=0.8508 精确值：精确值：cos (cos ( /6)=0.8660/6)=0.8660nhhhhxxxxxnhMxRnjjjjnnL324,2, 1:)(2011112)2)(1(4324)!1(1)(nnnnnnhhhhnxRLn1234)( xRn0 .30 .0 44 .7 1 0-34 .7 1 0-42021/6/16281.1.拉格朗日(Lagrange)多项式插值1.51.5 拉格朗日插值多项式的振荡?)(?)(x

17、RxLnnn55,11)(2xxxg63. 363. 3),()(limxxgxLnnTo MATLAB(runge)Runge现象现象: :-505-1.5-1-0.500.511.52y=1/(1+x2)n=2n=4n=6n=8n=102021/6/16292.2.分段线性插值分段线性插值xjxj-1xj+1x0 xn其它,0,)()()(1111110jjjjjjjjjjjnjjjnxxxxxxxxxxxxxxxlxlyxI计算量与n无关;n越大，误差越小.nnnxxxxgxI0),()(lim2021/6/16303. 3. 三次样条插值三次样条插值, 1,),()(1nixxxxsx

18、SiiiL,)()3), 1 ,0()()2), 1()()10223niiiiiiixxCxSniyxSnidxcxbxaxsLL) 1, 1()()(),()(),()(111 nixsxsxsxsxsxsiiiiiiiiiiii自然边界条件）(0)()()40 nxSxS)(,)4)3)2xSdcbaiiii)()(limxgxsn2021/6/16311. 1. 拉格朗日插值拉格朗日插值: :自编程序自编程序, ,如名为如名为 lagr1.m 的的M文件，文件， 第一行为第一行为 function y=lagr1(x0,y0,x) 输入输入: :节点节点x0,y0, 插值点插值点x (

19、 (均为数组，长度自定义均为数组，长度自定义) )）；）； 输出输出: :插值插值y ( (与与x同长度数组同长度数组) )）。）。 应用时输入应用时输入x0,y0,x后后, ,运行运行 y=lagr1(x0,y0,x)2. 2. 分段线性插值分段线性插值: :已有程序已有程序 y=interp1(x0,y0,x)3. 3. 三次样条插值三次样条插值: :已有程序已有程序 y=interp1(x0,y0,x, spline) 或或 y=spline(x0,y0,x)用MATLAB作插值计算2021/6/1632对表格给出的函数，求出没有给出的函数值 例例2表中是待加工零件下轮廓线的一组数据，现

20、需要得到x坐标每改变0.1时所对应的y的坐标. x03571112131415y01.21.72.12.01.81.21.01.6下面是关于插值的两条命令(专门用来解决这类问题)：y=interp1(x0,y0,x) 分段线性插值y=spline(x0,y0,x) 三次样条插值 2021/6/1633n其中x0,y0是已知的节点坐标，是同维向量。y对应于x处的插值。y与x是同维向量。n解决上述问题,我们可分两步:n一 用原始数据绘图作为选用插值方法的参考.n二 确定插值方法进行插值计算x0=0,3,5,7,9,11,12,13,14,15;y0=0,1.2,1.7,2.0,2.1,2.0,1.

21、8,1.2,1.0,1.6plot(x0,y0) %完成第一步工作x=0:0.1:15;y=interp1(x0,y0,x); %用分段线性插值完成第二步工作plot(x,y)y=spline(x0,y0,x); plot(x,y) %用三次样条插值完成第二步工作2021/6/1634n对y=1/(1+x2),-5x5，用n（=11）个节点（等分）作上述两种插值，用m（=21）个插值点（等分）作图，比较结果。n解：键入并运行如下命令nn=11;m=21;x=-5:10/(m-1):5;y=1./(1+x.2);nxo=-5:10/(n-1):5;yo=1./(1+xo.2);ny1=inter

22、p1(xo,yo,x);ny2=spline(xo,yo,x);nplot(x,y,r,x,y1,b,x,y2,k)练习：2021/6/1635例 3 一水库上游河段降暴雨.，根据预报测算上游流入水库的流量为Q(t) (102立方米/秒) ： t (时) 8 12 16 24 30 44 48 56 60Q( t ) 36 54 78 92 101 35 25 16 13 利用这个预报值估计14：30 和 20：30 时上游流入水库的流量。假设: 1 已知数据准确。 2 相邻两个时刻之间的流量没有突然的变化。2021/6/16362021/6/1637nt=8,12,16,24,30,44,4

23、8,56,60;nq=36,54,78,92,101,35,25,16,13;nt1=8:0.5:60;nq1=interp1(t,q,t1,linear);nplot(t,q,b,t1,q1);nhold on;nq2=interp1(t,q,t1,spline);nplot(t,q,b,t1,q1,t1,q2,r)nq1nq22021/6/16382021/6/1639拟合与插值的关系拟合与插值的关系 函数插值与曲线拟合都是要根据一组数据构造一个函数作函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者在数学方法上是完全不同为近似，由于近似的要求不同，二者在数

24、学方法上是完全不同的的 实例：实例：下面数据是某次实验所得，希望得到X和 f之间的关系？x124791 21 31 51 7f1 .53 .96 .611 .71 5 .61 8 .81 9 .62 0 .62 1 .1问题：问题：给定一批数据点，需确定满足特定要求的曲线或曲面解决方案：解决方案：若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，这就是数据拟合数据拟合，又称曲线拟合或曲面拟合若要求所求曲线（面）通过所给所有数据点，就是插值问题插值问题；2021/6/1640最临近插值、线性插值、样条插值与曲线拟合结果：最临近插值、线性插值、样条插值与曲线拟合结果：0246810121416180510152025已知数据点spline三次多项式插值0246810121416180510152025已知数据点linest三次多项式插值0246810121416180510152025已知数据点nearest三次多项式插值2021/6/1641即要求即要求 出二次多项式出二次多项式:3221)(axaxaxf中中 的的),(321aaaA 使得使得:最小 )(1112iiiyxf1、对