数学建模——数值计算方法总结学习教案.pptx

《数学建模——数值计算方法总结学习教案.pptx》由会员分享，可在线阅读，更多相关《数学建模——数值计算方法总结学习教案.pptx（87页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、会计学1数学数学(shxu)建模建模数值计算方法总结数值计算方法总结第一页，共87页。假设已获得某函数关系的成批离散实验数据或观测数据，假设已获得某函数关系的成批离散实验数据或观测数据，拟合问题就是为这样的大量离散数据建立对应的、近似的连拟合问题就是为这样的大量离散数据建立对应的、近似的连续模型的一种应用基础问题。所建立的模型的基本形式是一续模型的一种应用基础问题。所建立的模型的基本形式是一条曲线（一元条曲线（一元(y yun)曲线），称为拟合曲线或经验公式。曲线），称为拟合曲线或经验公式。通常采用通常采用“误差的平方和最小误差的平方和最小”的原则的原则(yunz)，即最，即最小二乘拟合问题。

2、小二乘拟合问题。它不要求目标模型（即拟合曲线它不要求目标模型（即拟合曲线(qxin)）精确地）精确地过已知的各离散点，只要求目标模型符合已知离散点分过已知的各离散点，只要求目标模型符合已知离散点分布的总体轮廓，并与已知离散点的误差按某种意义尽量布的总体轮廓，并与已知离散点的误差按某种意义尽量地小。地小。一、拟合问题一、拟合问题第2页/共87页第二页，共87页。拟拟 合合 问问 题题 引引 例例 1 1温度温度(wnd)t(0C)20.5 32.7 (wnd)t(0C)20.5 32.7 51.0 73.0 95.751.0 73.0 95.7电阻电阻R(R()765 826 873 )765

3、826 873 942 1032942 1032已知热敏电阻已知热敏电阻(r mn din z)数据：数据：求求600C600C时的电阻时的电阻(dinz)R(dinz)R。设设 R=at+ba,b为待定系数为待定系数第3页/共87页第三页，共87页。拟拟 合合 问问 题题 引引 例例 2 2 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)求血药浓度随时间求血药浓度随

4、时间(shjin)的变化规律的变化规律c(t).作半对数作半对数(du sh)坐标系坐标系(semilogy)下的图形下的图形第4页/共87页第四页，共87页。已知一组观测已知一组观测(gunc)数据：数据：要求在某特定函数类要求在某特定函数类 中寻找一个函数中寻找一个函数 作为作为的近似函数，使得二者在节点产生的残差的近似函数，使得二者在节点产生的残差按某种度量按某种度量(dling)标准为最标准为最小。小。常用常用(chn yn)原则：残差平方和最小原则：残差平方和最小曲曲 线线 拟拟 合合 问问 题题 的的 提提 法法第5页/共87页第五页，共87页。线性最小二乘拟合函数线性最小二乘拟合

5、函数(hnsh)的选取的选取 +=a1+a2x+=a1+a2x+a3x2+=a1+a2x+a3x2=a1+a2/x=aebx=ae-bx1.1.通过机理分析建立数学模型来确定通过机理分析建立数学模型来确定 ；2.2.将数据将数据 (xi,yi)i=1,n 作图，通过直观判断确定作图，通过直观判断确定 ：第6页/共87页第六页，共87页。曲线拟合问题最常用的解法曲线拟合问题最常用的解法(ji f)线性最小二乘法的基本线性最小二乘法的基本思路思路 第二步第二步:确定确定a1,a2,an 的准则的准则(zhnz)（最小二乘准则（最小二乘准则(zhnz)）：）：使使n个点（个点（xi,yi)与曲线与曲

6、线 y=(x)的距离的距离i 的平方和最小的平方和最小。记记问题问题(wnt)归结为，求归结为，求 a1,a2,an 使使 J(a1,a2,an)最小。最小。第一步:先选定一组函数先选定一组函数 使使其中其中 a1,a2,an 为待定系数。为待定系数。第7页/共87页第七页，共87页。方程组没有方程组没有(mi yu)通常意义下的解，这类方程组称为超定方通常意义下的解，这类方程组称为超定方当线性方程当线性方程(fngchng)组的方程组的方程(fngchng)个数多于未知数的个数时，个数多于未知数的个数时，设线性方程组为设线性方程组为程组或矛盾程组或矛盾(modn)方程组。方程组。最小二乘法的

7、求解：预备知识最小二乘法的求解：预备知识第8页/共87页第八页，共87页。若能找到一组向量若能找到一组向量(xingling)(xingling)令令最小，其中最小，其中(qzhng)(qzhng)使得使得(sh de)(sh de)则称则称 为该为该超定方程组的最小二乘解超定方程组的最小二乘解。由多元函数取极值的必要条件有由多元函数取极值的必要条件有求其最小值。求其最小值。第9页/共87页第九页，共87页。即即第10页/共87页第十页，共87页。故得故得第11页/共87页第十一页，共87页。即即称为称为(chn wi)正则方程组。正则方程组。该方程组的解即为超定方程组的最小二乘解。该方程组的

8、解即为超定方程组的最小二乘解。第12页/共87页第十二页，共87页。(2)求解求解(qi ji)正则方程组得最小二乘解。正则方程组得最小二乘解。用最小二乘法用最小二乘法(chngf)解超定方程组的步骤：解超定方程组的步骤：(1)计算计算 和和 ，得正则方程组，得正则方程组 。第13页/共87页第十三页，共87页。例例解得最小二乘解为解得最小二乘解为得方程组得方程组解：解：已知试验数据，用最小二乘法求拟合已知试验数据，用最小二乘法求拟合(n h)直线直线0.00.20.40.60.80.91.92.83.34.2故得拟合故得拟合(n h)直线直线第14页/共87页第十四页，共87页。可线性化模型

9、可线性化模型(mxng)的最小二乘的最小二乘拟合拟合很多实际很多实际(shj)问题中，变量间并非线性关系，但拟合问题中，变量间并非线性关系，但拟合曲线可视为曲线可视为 的形式，的形式，指数函数指数函数(zh sh hn sh)如双曲线如双曲线即将非线性化问题转化为线性问题。即将非线性化问题转化为线性问题。令令则有则有第15页/共87页第十五页，共87页。例例给定如下观测数据，试用指数曲线给定如下观测数据，试用指数曲线 进行拟合。进行拟合。1.01.251.51.752.05.15.796.537.458.46解：解：令令则有则有1.01.251.51.752.01.629 1.756 1.87

10、6 2.008 2.135故故解此超定方程组得解此超定方程组得则拟合则拟合(n h)曲线为曲线为第16页/共87页第十六页，共87页。多变量多变量(binling)数数据拟合据拟合有时变量间关系为多元有时变量间关系为多元(du yun)函数关系，函数关系，有如下有如下(rxi)观测数据观测数据观测次数12m假定变量假定变量y与与n个变量个变量xi间为线性关系，间为线性关系，可设拟合方程为可设拟合方程为第17页/共87页第十七页，共87页。第第i组观测组观测(gunc)数据对应的残差为数据对应的残差为下面考虑用最小二乘原理下面考虑用最小二乘原理(yunl)确定拟合方程的系数确定拟合方程的系数ai

11、。按照最小二乘原理，应使按照最小二乘原理，应使 最小。最小。令令求解求解(qi ji)该超定方程组的该超定方程组的最小二乘解即可。最小二乘解即可。第18页/共87页第十八页，共87页。用用MATLAB 解拟合解拟合(n h)问题问题1 1、线性最小二乘拟合、线性最小二乘拟合(n h)(n h)2 2、非线性最小二乘拟合、非线性最小二乘拟合(n(n h)h)第19页/共87页第十九页，共87页。用用MATLABMATLAB作线性最小二乘拟合作线性最小二乘拟合(n h)(n h)1.1.作多项式作多项式f(x)=a1xm+amx+am+1f(x)=a1xm+amx+am+1拟合拟合(n h),(n

12、 h),可利用可利用已有命令已有命令:a=polyfit(x,y,m)3.3.对超定方程组对超定方程组可得最小二乘意义下的解。可得最小二乘意义下的解。，用用2.2.多项式在多项式在x x处的值处的值y y的计算的计算(j sun)(j sun)命令：命令：y=polyvaly=polyval（a a，x x）输出拟合多项式系数输出拟合多项式系数a=a1,am,am+1（数组）数组）输入同长度输入同长度数组数组X，Y拟合多项式拟合多项式次数次数第20页/共87页第二十页，共87页。即要求即要求 出二次多项式出二次多项式:中中 的的使得使得:例例 对下面一组数据对下面一组数据(shj)作二次多项式

13、拟合作二次多项式拟合第21页/共87页第二十一页，共87页。1）输入）输入(shr)命令：命令：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=Ry解法解法1 1解超定方程解超定方程(fngchng)(fngchng)的方法的方法2）计算结果）计算结果：=-9.8108,20.1293,-0.0317第22页/共87页第二十二页，共87页。第23页/共87页第二十三页，共87页。1）输入命令：）输入命令：x=0:0.1:1;y=-0.447,1.978,3.

14、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)%作出数据点和拟合作出数据点和拟合(n h)曲线的图形曲线的图形2）计算结果：）计算结果：=-9.8108,20.1293,-0.0317解法解法2用多项式拟合用多项式拟合(n h)的命令的命令第24页/共87页第二十四页，共87页。1.lsqcurvefit1.lsqcurvefit已知数据已知数据(shj)(shj)点：点：xdata=xdata=（xdata1xdata1，xdata2xdata2，xdat

15、anxdatan）ydata=ydata=（ydata1ydata1，ydata2ydata2，ydatanydatan）用用MATLABMATLAB作非线性最小二乘拟合作非线性最小二乘拟合(n h)(n h)两个求非线性最小二乘拟合两个求非线性最小二乘拟合(n h)(n h)的函数：的函数：lsqcurvefitlsqcurvefit、lsqnonlinlsqnonlin。相同点和不同点：两个命令都要先建立相同点和不同点：两个命令都要先建立M-M-文件文件fun.mfun.m，定义函数，定义函数f(x)f(x)，但定义，但定义f(x)f(x)的方式不同。的方式不同。lsqcurvefitls

16、qcurvefit用以求含参量用以求含参量x x（向量）的向量值函数（向量）的向量值函数F(x,xdata)=F(x,xdata)=（F F（x x，xdataxdata1 1），），F F（x x，xdataxdatan n）T T使得使得 第25页/共87页第二十五页，共87页。输入输入(shr)(shr)格式格式:(1)x=lsqcurvefit(fun,x0,xdata,ydata);(1)x=lsqcurvefit(fun,x0,xdata,ydata);(2)x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub);(2)x=lsqcurvefit(fun,x

17、0,xdata,ydata,lb,ub);(3)x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,(3)x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options);options);(4)x,resnorm=lsqcurvefit (4)x,resnorm=lsqcurvefit(fun,x0,xdata,ydata,);(fun,x0,xdata,ydata,);(5)x,resnorm,residual=lsqcurvefit (5)x,resnorm,residual=lsqcurvefit(fun,x0,xdata,yda

18、ta,);(fun,x0,xdata,ydata,);fun是一个事先建立的是一个事先建立的定义函数定义函数F(x,xdata)的的M-文件文件,自变量为自变量为x和和xdata说明说明(shumng)(shumng)：x=lsqcurvefit(fun,x0,xdata,ydata,options);x=lsqcurvefit(fun,x0,xdata,ydata,options);迭代初值迭代初值已知数据点已知数据点选项见无选项见无约束优化约束优化第26页/共87页第二十六页，共87页。lsqnonlin用以求含参量用以求含参量x x（向量）的向量值函数（向量）的向量值函数 f(x)f(x

19、)=(f=(f1 1(x),f(x),f2 2(x),(x),f,fn n(x)(x)T T ，使得，使得 最小。最小。其中其中 fi（x）=f（x，xdatai，ydatai）=F(x,xdatai)-ydatai2.lsqnonlin已知数据已知数据(shj)(shj)点：点：xdata=xdata=（xdata1xdata1，xdata2xdata2，xdatanxdatan）ydata=ydata=（ydata1ydata1，ydata2ydata2，ydatanydatan）第27页/共87页第二十七页，共87页。输入输入(shr)格式：格式：1）x=lsqnonlin（fun，x0

20、）；）；2）x=lsqnonlin（fun，x0，lb，ub）；）；3）x=lsqnonlin（fun，x0，lb，ub，options）；）；4）x，resnorm=lsqnonlin（fun，x0，）；）；5）x，resnorm，residual=lsqnonlin（fun，x0，）；）；说明说明(shumng)：x=lsqnonlin（fun，x0，options）；）；fun是一个事先建立的是一个事先建立的定义函数定义函数f(x)的的M-文件，文件，自变量为自变量为x迭代初值迭代初值选项见无选项见无约束优化约束优化第28页/共87页第二十八页，共87页。例例2 用下面一组数据用下面一组

21、数据(shj)拟合拟合 中的参数中的参数a，b，k该问题该问题(wnt)即解的最优化问题即解的最优化问题(wnt)：第29页/共87页第二十九页，共87页。1）编写(binxi)M-文件 curvefun1.m function f=curvefun1(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata)%其中 x(1)=a;x(2)=b；x(3)=k;2）输入）输入(shr)命令命令tdata=100:100:1000cdata=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;x0=0.2,0.05

22、,0.05;x=lsqcurvefit(curvefun1,x0,tdata,cdata)f=curvefun1(x,tdata)F(x，tdata)=，x=(a，b，k)解法解法(ji f)1.(ji f)1.用命令用命令lsqcurvefitlsqcurvefit第30页/共87页第三十页，共87页。3 3）运算）运算(yn sun)(yn sun)结果：结果：f=0.0043 0.0051 0.0056 0.0059 0.0061 f=0.0043 0.0051 0.0056 0.0059 0.0061 0.0062 0.0062 0.0063 0.0063 0.0063 0.0062

23、0.0062 0.0063 0.0063 0.0063 x=0.0063 -0.0034 0.2542 x=0.0063 -0.0034 0.25424）结论）结论(jiln)：a=0.0063,b=-0.0034,k=0.2542第31页/共87页第三十一页，共87页。1）编写(binxi)M-文件 curvefun2.m function f=curvefun2(x)tdata=100:100:1000;cdata=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;f=x(1)+x(2)*exp(-0.02*x(3)*tdata

24、)-cdata2）输入）输入(shr)命令命令:x0=0.2,0.05,0.05;x=lsqnonlin(curvefun2,x0)f=curvefun2(x)函数函数(hnsh)curvefun2的自变量的自变量是是x，cdata和和tdata是已知参数，是已知参数，故应将故应将cdata tdata的值写在的值写在curvefun2.m中中解法解法 2 2 用命令用命令lsqnonlinlsqnonlin x=x=（a a，b b，k k）第32页/共87页第三十二页，共87页。3 3）运算）运算(yn sun)(yn sun)结果为结果为 f=1.0e-003*f=1.0e-003*（0

25、.2322 -0.1243 -0.2495 -0.2322 -0.1243 -0.2495 -0.2413 0.2413-0.1668 -0.0724 0.0241 0.1159 0.2030 -0.1668 -0.0724 0.0241 0.1159 0.2030 0.27920.2792）x=0.0063 -0.0034 0.2542 x=0.0063 -0.0034 0.2542可以看出可以看出(kn ch)，两个命令的计算结果是相同的。，两个命令的计算结果是相同的。4）结论）结论(jiln)：即拟合得：即拟合得a=0.0063 b=-0.0034 k=0.2542第33页/共87页第三

26、十三页，共87页。设有一个年产设有一个年产240吨的食品加工厂吨的食品加工厂,需要统计其生产需要统计其生产(shngchn)费费用用,但由于该厂各项资料不全但由于该厂各项资料不全,无法统计。在这种情况下无法统计。在这种情况下,统计部门统计部门收集了设备、生产收集了设备、生产(shngchn)能力和该厂大致相同的五个食品加工能力和该厂大致相同的五个食品加工厂的产量与生产厂的产量与生产(shngchn)费用资料如下表：费用资料如下表：如何根据已知五个食品加工厂的资料较准确地估计该厂的生产如何根据已知五个食品加工厂的资料较准确地估计该厂的生产(shngchn)费用呢？费用呢？引例引例(yn l)插值

27、问题插值问题(wnt)第34页/共87页第三十四页，共87页。第35页/共87页第三十五页，共87页。（2）有的函数）有的函数(hnsh)甚至没有表达式，只是一种表格函数甚至没有表达式，只是一种表格函数(hnsh)，（1）如果）如果(rgu)函数表达式本身比较复杂，且需要多次函数表达式本身比较复杂，且需要多次实际实际(shj)问题中经常要涉及到函数值的计算问题：问题中经常要涉及到函数值的计算问题：重复计算时，计算量会很大；重复计算时，计算量会很大；方便且表达简单的函数来近似代替，这就是方便且表达简单的函数来近似代替，这就是插值插值问题。问题。实际对于这两种情况，我们都需要寻找一个计算实际对于这

28、两种情况，我们都需要寻找一个计算而我们需要的函数值可能不在该表格中。而我们需要的函数值可能不在该表格中。问题问题：插值与拟合的区别？：插值与拟合的区别？第36页/共87页第三十六页，共87页。一一 维维 插插 值值一、插值的定义一、插值的定义(dngy)二、插值的方法二、插值的方法(fngf)三、用三、用Matlab解插值问题解插值问题(wnt)拉格朗日插值拉格朗日插值牛顿插值牛顿插值反插值反插值三次样条插值三次样条插值分段插值法分段插值法第37页/共87页第三十七页，共87页。二维插值二维插值一、二维插值定义一、二维插值定义(dngy)(dngy)二、网格二、网格(wn)(wn)节点插值法节

29、点插值法三、用三、用MatlabMatlab解插值问题解插值问题(wnt)(wnt)最邻近插值最邻近插值分片线性插值分片线性插值双线性插值双线性插值网格节点数据的插值网格节点数据的插值散点数据的插值散点数据的插值第38页/共87页第三十八页，共87页。一维插值的定义一维插值的定义(dngy)(dngy)已知已知 n+1个节点个节点其中其中互不相同，不妨设互不相同，不妨设求任一插值点求任一插值点处的插值处的插值第39页/共87页第三十九页，共87页。构造一个构造一个(相对简单的相对简单的)函数函数通过全部节点通过全部节点,即即再用再用计算插值，计算插值，即即第40页/共87页第四十页，共87页。

30、代数代数(dish)插值的唯一性插值的唯一性是惟一是惟一(wiy)的。的。定理定理个不同节点，满足插值条件个不同节点，满足插值条件的的n次插值多项式次插值多项式当选择代数多项式作为当选择代数多项式作为(zuwi)插值函数时，称为代数多插值函数时，称为代数多项式插值。项式插值。代数插值代数插值第41页/共87页第四十一页，共87页。一、线性插值一、线性插值一、线性插值一、线性插值（n n=1=1，一次插值），一次插值），一次插值），一次插值）求解求解(qi ji)L1(x)=a1 x+a0已知已知使得使得(sh de)L1(xi)=yi.(i=0,1)如果令如果令则称则称 l0(x),l1(x)

31、为为x0,x1上的线性插值基函数上的线性插值基函数(hnsh)。这时这时根据点斜式得到根据点斜式得到并称其为一次并称其为一次Lagrange插值多项式。插值多项式。f(x)L1(x)=y0l0(x)+y1 l1(x)第42页/共87页第四十二页，共87页。二、抛物线插值二、抛物线插值二、抛物线插值二、抛物线插值（n n=2=2，二次插值），二次插值），二次插值），二次插值）求解求解(qi ji)L2(x)=a2x2+a1 x+a0使得使得(sh de)L2(xi)=yi ,i=0,1,2.已知已知第43页/共87页第四十三页，共87页。抛物插值抛物插值仿照线性函数仿照线性函数(hnsh)的构造

32、方法，构造的构造方法，构造且且其中要求其中要求 均为二次多项式。均为二次多项式。设设即即由由求求第44页/共87页第四十四页，共87页。故故同理同理第45页/共87页第四十五页，共87页。例题例题(lt)求线性插值函数求线性插值函数(hnsh)。已知已知 的函数的函数(hnsh)表表 解：解：第46页/共87页第四十六页，共87页。多项式可使用类似多项式可使用类似(li s)方法。方法。分析分析(fnx)由由 个不同节点个不同节点 构造构造 次次其中其中(qzhng)且且上述多项式称为上述多项式称为拉格朗日插值多项式拉格朗日插值多项式。拉格朗日插值多项式拉格朗日插值多项式第47页/共87页第四

33、十七页，共87页。插值余项和误差插值余项和误差(wch)估计估计且依赖于且依赖于 。其中其中设设f(x)在在a,b有有n+1阶导数，则阶导数，则n次插值多次插值多定理定理项式项式 对任意对任意 ，有插值余项，有插值余项第48页/共87页第四十八页，共87页。拉格朗日多项式插值的这种振荡(zhndng)现象叫 Runge现象 采用拉格朗日多项式插值：选取(xunq)不同插值节点个数n+1，其中n为插值多项式的次数，当n分别取2,4,6,8,10时，绘出插值结果图形.例例第49页/共87页第四十九页，共87页。第50页/共87页第五十页，共87页。分段分段(fn(fn dun)dun)线性插线性插

34、值值xjxj-1xj+1x0 xnxoy第51页/共87页第五十一页，共87页。例例用分段用分段(fn dun)线性插值法求插值线性插值法求插值,并观察插值误差并观察插值误差.在在-6,6中平均中平均(pngjn)选取选取41个点作插值个点作插值,结结果如图示果如图示第52页/共87页第五十二页，共87页。比分比分(b fn)段线性插值更光滑。段线性插值更光滑。xyxi-1 xiab 在数学上，光滑程度的定量描述是：函数(曲线)的k阶导数存在且连续，则称该曲线具有k阶光滑性。光滑性的阶次越高，则越光滑。是否(sh fu)存在较低次的分段多项式达到较高阶光滑性的方法？三次样条插值就是一个很好的例

35、子。三次三次(sn c)(sn c)样条插值样条插值第53页/共87页第五十三页，共87页。所谓样条函数，从数学上说，就是按照一定光滑性所谓样条函数，从数学上说，就是按照一定光滑性它在每个内节点上具有直到它在每个内节点上具有直到m-1m-1阶连续导数。阶连续导数。要求要求“装配装配”起来的分段多项式。起来的分段多项式。它在每个子区间内是它在每个子区间内是m m次多项式，次多项式，三次三次(sn c)(sn c)样条插值样条插值三次三次(sn c)样条函数满足：样条函数满足：在每个子区间在每个子区间 上是一个三次多项式。上是一个三次多项式。(1)(2)在每个内节点上具有直到二阶的连续导数。在每个

36、内节点上具有直到二阶的连续导数。(3)上可设上可设故在每一个小区间故在每一个小区间求解求解(qi ji)使用三弯矩法使用三弯矩法第54页/共87页第五十四页，共87页。例例用三次用三次(sn c)样条插值选取样条插值选取11个基点计个基点计算插值算插值第55页/共87页第五十五页，共87页。用用MATLABMATLAB作插值计算作插值计算(j sun)(j sun)一维插值函数一维插值函数(hnsh)(hnsh)：yi=interp1(x，y，xi，method)插值方法插值方法被插值点被插值点插值节点插值节点xixi处的插处的插值结果值结果nearest：最邻近插值：最邻近插值linear：

37、线性插值；线性插值；spline：三次样条插值；三次样条插值；cubic：立方立方(lfng)插值。插值。缺省时：缺省时：分段线性插值。分段线性插值。注意：所有的插值方法都要求注意：所有的插值方法都要求x是单调的，并且是单调的，并且xi不能够超不能够超过过x的范围。的范围。第56页/共87页第五十六页，共87页。例：在例：在1-121-12的的1111小时小时(xiosh)(xiosh)内，每隔内，每隔1 1小时小时(xiosh)(xiosh)测量一次温度，测得的温度依次为：测量一次温度，测得的温度依次为：5 5，8 8，9 9，1515，2525，2929，3131，3030，2222，25

38、25，2727，2424。试估计每。试估计每隔隔1/101/10小时小时(xiosh)(xiosh)的温度值。的温度值。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);(直接(zhji)输出数据将是很多的)plot(hours,temps,+,h,t,hours,temps,r:)%作图xlabel(Hour),ylabel(Degrees Celsius)第57页/共87页第五十七页，共87页。第58页/共87页第五十八页，共87页。xy机翼下轮廓线例例

39、已知飞机已知飞机(fij)(fij)下轮廓线上数据如下，求下轮廓线上数据如下，求x x每改变每改变0.10.1时的时的y y值。值。第59页/共87页第五十九页，共87页。第60页/共87页第六十页，共87页。二维插值的定义二维插值的定义(dngy)xyO O第一种（网格第一种（网格(wn)节点）：节点）：第61页/共87页第六十一页，共87页。已知已知 mn个节点个节点(ji din)其中其中互不相同，不妨设互不相同，不妨设 构造一个二元函数构造一个二元函数通过全部已知节点通过全部已知节点,即即再用再用计算插值，计算插值，即即第62页/共87页第六十二页，共87页。第二种（散乱第二种（散乱(

40、sn lun)节节点）：点）：yx0 0第63页/共87页第六十三页，共87页。已知已知n个节点个节点其中其中互不相同，互不相同，构造一个二元函数构造一个二元函数通过全部已知节点通过全部已知节点,即即再用再用计算插值，计算插值，即即第64页/共87页第六十四页，共87页。注意：最邻近插值一般不连续。具有注意：最邻近插值一般不连续。具有(jyu)连续性连续性的最简单的插值是分片线性插值。的最简单的插值是分片线性插值。最邻近最邻近(ln(ln jn)jn)插值插值x y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O O 二维或高维情形的最邻近(ln jn)插值，与被插值点最邻近(ln

41、jn)的节点的函数值即为所求。第65页/共87页第六十五页，共87页。将 四 个 插 值 点（矩 形 的 四 个 顶 点）处 的 函 数(hnsh)值依次简记为：分片分片(fn pin)(fn pin)线线性插值性插值xy (xi,yj)(xi,yj+1)(xi+1,yj)(xi+1,yj+1)O Of(xi,yj)=f1，f(xi+1,yj)=f2，f(xi+1,yj+1)=f3，f(xi,yj+1)=f4第66页/共87页第六十六页，共87页。插值函数(hnsh)为：第二片(上三角形区域(qy)：(x,y)满足插值函数(hnsh)为：注意注意：(x,y)当然应该是在插值节点所形成的矩形区域

42、内。显然，分片线性插值函数是连续的；分两片的函数表达式如下：第一片(下三角形区域)：(x,y)满足第67页/共87页第六十七页，共87页。双线性插值是一片(y pin)一片(y pin)的空间二次曲面构成。双线性插值函数的形式如下：其中(qzhng)有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。双线性插值双线性插值x y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O O第68页/共87页第六十八页，共87页。要求要求x0,y0 x0,y0单调单调(dndio)(dndio)；x x，y y可取为矩可取为矩阵，或阵，或x x取行向

43、量，取行向量，y y取为列向量，取为列向量，x,yx,y的值分的值分别不能超出别不能超出x0,y0 x0,y0的范围。的范围。z=interp2(x0,y0,z0,x,y,method)被插值点插值方法用用MATLABMATLAB作网格作网格(wn)(wn)节点数据节点数据的插值的插值插值节点被插值点的函数值nearest 最邻近最邻近(ln jn)插值插值linear 双线性插值双线性插值cubic 双三次插值双三次插值缺省时缺省时,双线性插值双线性插值第69页/共87页第六十九页，共87页。例：测得平板表面例：测得平板表面3*53*5网格点处的温度分别为：网格点处的温度分别为：82 81

44、80 82 84 82 81 80 82 84 79 63 61 65 81 79 63 61 65 81 84 84 82 85 86 84 84 82 85 86 试作出平板表面的温度分布试作出平板表面的温度分布(fnb)(fnb)曲面曲面z=f(x,y)z=f(x,y)的图形。的图形。输入以下(yxi)命令：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)1.先在三维坐标画出原始数据，画出粗糙的温度(wnd)分布曲图.第70页/共87页第七十页，共87页。第71页/共87页第七十一

45、页，共87页。再输入以下命令(mng lng):xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi,yi,cubic);mesh(xi,yi,zi)画出插值后的温度分布曲面图.2以平滑数据,在x、y方向(fngxing)上每隔0.2个单位的地方进行插值.第72页/共87页第七十二页，共87页。第73页/共87页第七十三页，共87页。通过此例对最近邻点插值、双线性插值方法和双三次(sn c)插值方法的插值效果进行比较。第74页/共87页第七十四页，共87页。原始数据图原始数据图第75页/共87页第七十五页，共87页。最邻近最邻近(ln(ln jn)jn)插

46、值插值第76页/共87页第七十六页，共87页。双线性插值第77页/共87页第七十七页，共87页。双三次(sn c)插值第78页/共87页第七十八页，共87页。插值函数(hnsh)griddata格式为:cz=griddata（x，y，z，cx，cy，method）用用MATLABMATLAB作散点数据作散点数据(shj)(shj)的插值计的插值计算算 要求要求(yoqi)cx(yoqi)cx取行向量，取行向量，cycy取为列向量。取为列向量。被插值点插值方法插值节点被插值点的函数值nearestnearest 最邻近插值最邻近插值linearlinear 双线性插值双线性插值cubiccubi

47、c 双三次插值双三次插值v4-Matlab提供的插值方法提供的插值方法缺省时缺省时,双线性插值双线性插值第79页/共87页第七十九页，共87页。一室模型：将整个机体看作一个房室，称中心室，一室模型：将整个机体看作一个房室，称中心室，室内血药浓度是均匀室内血药浓度是均匀(jnyn)的。快速静脉注射后，的。快速静脉注射后，浓度立即上升；然后迅速下降。当浓度太低时，达浓度立即上升；然后迅速下降。当浓度太低时，达不到预期的治疗效果；当浓度太高，又可能导致药不到预期的治疗效果；当浓度太高，又可能导致药物中毒或副作用太强。临床上，每种药物有一个最物中毒或副作用太强。临床上，每种药物有一个最小有效浓度小有效

48、浓度c1和一个最大有效浓度和一个最大有效浓度c2。设计给药方。设计给药方案时，要使血药浓度案时，要使血药浓度 保持在保持在c1c2之间。本题设之间。本题设c1=10，c2=25(ug/ml).拟拟 合合 问问 题题 实实 例例给药方案给药方案 一种新药用于临床之前一种新药用于临床之前(zhqin)，必须设计给药方案，必须设计给药方案.药物进入机体后血液输送到全身，在这个过程中不药物进入机体后血液输送到全身，在这个过程中不断地被吸收断地被吸收(xshu)、分布、代谢，最终排出体外，药物、分布、代谢，最终排出体外，药物在血液中的浓度，即单位体积血液中的药物含量，称为在血液中的浓度，即单位体积血液中

49、的药物含量，称为血药浓度。血药浓度。第80页/共87页第八十页，共87页。在实验方面在实验方面,对某人用快速静脉注射方式一次对某人用快速静脉注射方式一次注入该药物注入该药物300mg后后,在一定时刻在一定时刻t(小时小时)采集血药采集血药,测得血药浓度测得血药浓度c(ug/ml)如下表如下表: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 要设计给药方案要设计给药方案,必须知道给药后血药浓度随时间必须知道给药后血药浓度随时间变化变化(binhu)的规律。从实验和理论两方面着

50、手：的规律。从实验和理论两方面着手：第81页/共87页第八十一页，共87页。给药方案给药方案(fng n)1.在快速静脉注射的给药方式下，研究血药浓度在快速静脉注射的给药方式下，研究血药浓度（单位体积（单位体积(tj)血液中的药物含量）的变化规律。血液中的药物含量）的变化规律。tc2cc10问问题题(wnt)2.给定药物的最小有效浓度和最大治疗浓度，设计给药方给定药物的最小有效浓度和最大治疗浓度，设计给药方案：每次注射剂量多大；间隔时间多长。案：每次注射剂量多大；间隔时间多长。分分析析 理论：用一室模型研究血理论：用一室模型研究血药浓度变化规律药浓度变化规律 实验：对血药浓实验：对血药浓度数据