数学建模课件——数值计算方法

1 6/25/2024 数学建模教程数学建模教程拟拟 合与合与 插插 值值 2 6/25/2024 q在大量的应用领域中，人们经常面临这样的问题：在大量的应用领域中，人们经常面临这样的问题：给定一批数据点，需确定满足特定要求的曲线或给定一批数据点，需确定满足特定要求的曲线或曲面。对这个问题有两种方法。曲面。对这个问题有两种方法。q一种是一种是插值法插值法，数据假定是正确的，要求以某种方法描述数，数据假定是正确的，要求以某种方法描述数据点之间所发生的情况。据点之间所发生的情况。q另一种方法是另一种方法是曲线拟合或回归曲线拟合或回归。人们设法找出某条光滑曲线，。人们设法找出某条光滑曲线，它最佳地拟合数据，但不必要经过任何数据点。它最佳地拟合数据，但不必要经过任何数据点。q 本专题的目的之一是：了解插值和拟合的基本内容及本专题的目的之一是：了解插值和拟合的基本内容及方法；方法；函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。由于近似的要求不同，二者的数学方法上是完全不同的。3 6/25/2024 假设已获得某函数关系的假设已获得某函数关系的成批离散实验数据成批离散实验数据或或观观测数据测数据，拟合问题就是为这样的大量离散数据建立对，拟合问题就是为这样的大量离散数据建立对应的、近似的连续模型的一种应用基础问题。所建立应的、近似的连续模型的一种应用基础问题。所建立的模型的基本形式是一条曲线（一元曲线），称为的模型的基本形式是一条曲线（一元曲线），称为拟拟合曲线或经验公式。合曲线或经验公式。通常采用通常采用“误差的平方和最小误差的平方和最小”的原则，即的原则，即最小最小二乘拟合问题。二乘拟合问题。它它不要求不要求目标模型（即拟合曲线）目标模型（即拟合曲线）精确地精确地过已知过已知的各离散点，只要求目标模型的各离散点，只要求目标模型符合符合已知离散点分布的已知离散点分布的总体轮廓总体轮廓，并与已知离散点的误差，并与已知离散点的误差按某种意义按某种意义尽量地尽量地小小。一、拟合问题一、拟合问题4 6/25/2024 拟拟 合合 问问 题题 引引 例例 1 1温度温度t(0C)20.5 32.7 51.0 73.0 95.7电阻电阻R()765 826 873 942 1032已知热敏电阻数据：已知热敏电阻数据：求求60600C时的电阻时的电阻R。设设 R=at+ba,b为待定系数为待定系数5 6/25/2024 拟拟 合合 问问 题题 引引 例例 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)求血药浓度随时间的变化规律求血药浓度随时间的变化规律c(t).作半对数坐标系作半对数坐标系(semilogy)下的图形下的图形6 6/25/2024 已知一组观测数据：已知一组观测数据：要求在某特定函数类要求在某特定函数类 中寻找一个函数中寻找一个函数 作为作为的近似函数，使得二者在节点产生的残差的近似函数，使得二者在节点产生的残差按按某种度量标准为最小某种度量标准为最小。常用原则常用原则：残差平方和最小：残差平方和最小曲曲 线线 拟拟 合合 问问 题题 的的 提提 法法7 6/25/2024 线性最小二乘拟合函数线性最小二乘拟合函数的选取的选取 +=a1+a2x+=a1+a2x+a3x2+=a1+a2x+a3x2=a1+a2/x=aebx=ae-bx1.1.通过机理分析建立数学模型来确定通过机理分析建立数学模型来确定 ；2.2.将数据将数据 (xi,yi)i=1,n 作图，通过直观判断确定作图，通过直观判断确定 ：8 6/25/2024 曲线拟合问题最常用的解法曲线拟合问题最常用的解法线性最小二乘法的基本思路线性最小二乘法的基本思路 第二步:确定确定a1,a2,an 的的准则（最小二乘准则）：准则（最小二乘准则）：使使n个点个点（xi,yi)与与曲线曲线 y=(x)的的距离距离 i 的平方和最小的平方和最小。记记问题归结为，求问题归结为，求 a1,a2,an 使使 J(a1,a2,an)最小。最小。第一步:先选定一组函数先选定一组函数 使使其中其中 a1,a2,an 为待定系数。为待定系数。9 6/25/2024 方程组没有通常意义下的解，这类方程组称为方程组没有通常意义下的解，这类方程组称为超定方超定方当线性方程组的方程个数多于未知数的个数时，当线性方程组的方程个数多于未知数的个数时，设线性方程组为设线性方程组为程组程组或或矛盾方程组矛盾方程组。最小二乘法的求解：预备知识最小二乘法的求解：预备知识10 6/25/2024 若能找到一组向量若能找到一组向量令令最小，其中最小，其中使得使得则称则称 为该为该超定方程组的最小二乘解超定方程组的最小二乘解。由多元函数取极值的必要条件有由多元函数取极值的必要条件有求其最小值。求其最小值。11 6/25/2024 即即12 6/25/2024 故得故得13 6/25/2024 即即称为正则方程组。称为正则方程组。该方程组的解即为超定方程组的最小二乘解。该方程组的解即为超定方程组的最小二乘解。14 6/25/2024 (2)求解正求解正则方程方程组得最小二乘解。得最小二乘解。用最小二乘法解超定方程用最小二乘法解超定方程组的步的步骤：(1)计算算 和和 ，得正，得正则方程方程组 。15 6/25/2024 例例解得最小二乘解解得最小二乘解为得方程得方程组解：解：已知已知试验数据，用最小二乘法求数据，用最小二乘法求拟合直合直线0.00.20.40.60.80.91.92.83.34.2故得故得拟合直合直线16 6/25/2024 可线性化模型的最小二乘拟合可线性化模型的最小二乘拟合很多很多实际问题中，中，变量量间并非并非线性关系，但性关系，但拟合合曲曲线可可视为 的形式，的形式，指数函数指数函数如双曲如双曲线即将非即将非线性化性化问题转化化为线性性问题。令令则有有17 6/25/2024 例例给定如下观测数据，试用指数曲线给定如下观测数据，试用指数曲线 进行拟合。进行拟合。1.01.251.51.752.05.15.796.537.458.46解解：令令则有则有1.01.251.51.752.01.6291.7561.8762.0082.135故故解此超定方程组得解此超定方程组得则拟合曲线为则拟合曲线为18 6/25/2024 多变量数据拟合多变量数据拟合有时变量间关系为多元函数关系，有时变量间关系为多元函数关系，有如下观测数据有如下观测数据观测次数12m假定变量假定变量y与与n个变量个变量xi间为线性关系，间为线性关系，可设拟合方程为可设拟合方程为19 6/25/2024 第第i组观测数据对应的残差为组观测数据对应的残差为下面考虑用最小二乘原理确定拟合方程的系数下面考虑用最小二乘原理确定拟合方程的系数ai。按照最小二乘原理，应使按照最小二乘原理，应使 最小。最小。令令求解该超定方程组的求解该超定方程组的最小二乘解即可。最小二乘解即可。20 6/25/2024 用用MATLAB解拟合问题解拟合问题1 1、线性最小二乘拟合、线性最小二乘拟合2 2、非线性最小二乘拟合、非线性最小二乘拟合21 6/25/2024 用用MATLAB作线性最小二乘拟合作线性最小二乘拟合1.1.作作多多项式项式f(x)=a1xm+amx+am+1拟合拟合,可利用已有命令可利用已有命令:a=polyfit(x,y,m)3.3.对超定方程组对超定方程组可得最小二乘意义下的解。可得最小二乘意义下的解。，用，用2.2.多项式在多项式在x x处的值处的值y y的计算命令的计算命令：y=y=polyvalpolyval（a a，x x）输出拟合多项式系数输出拟合多项式系数a=a1,am,am+1（数组）数组）输入同长度输入同长度数组数组X，Y拟合多项式拟合多项式次数次数22 6/25/2024 即要求即要求 出二次多项式出二次多项式:中中 的的使得使得:例例 对下面一组数据作二次多项式拟合对下面一组数据作二次多项式拟合23 6/25/2024 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=Ry解法解法1 1解超定方程的方法解超定方程的方法2）计算结果）计算结果：=-9.8108,20.1293,-0.031724 6/25/2024 25 6/25/2024 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）计算结果：）计算结果：=-9.8108,20.1293,-0.0317解法解法2用多项式拟合的命令用多项式拟合的命令26 6/25/2024 1.1.lsqcurvefitlsqcurvefit已知数据点数据点：xdataxdata=（xdata1，xdata2，xdataxdatan n）ydataydata=（ydataydata1 1，ydataydata2 2，ydataydatan n）用用MATLAB作非线性最小二乘拟合作非线性最小二乘拟合两个求非线性最小二乘拟合的函数：两个求非线性最小二乘拟合的函数：lsqcurvefit、lsqnonlin。相同点和不同点：两相同点和不同点：两个命令都要先建立个命令都要先建立M-M-文件文件fun.mfun.m，定义函，定义函数数f(x)f(x)，但定义但定义f(x)f(x)的方式不同的方式不同。lsqcurvefitlsqcurvefit用以求含参量用以求含参量x x（向量）的向量值函数向量）的向量值函数F(x,xdataF(x,xdata)=)=（F F（x x，xdataxdata1 1），），F F（x x，xdataxdatan n）T T使得使得 27 6/25/2024 输入格式输入格式:(1)x=lsqcurvefit(fun,x0,xdata,ydata);(2)x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub);(3)x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options);(4)x,resnorm=lsqcurvefit(fun,x0,xdata,ydata,);(5)x,resnorm,residual=lsqcurvefit(fun,x0,xdata,ydata,);fun是一个事先建立的是一个事先建立的定义函数定义函数F(x,xdata)的的M-文件文件,自变量为自变量为x和和xdata说明：说明：x=lsqcurvefit(fun,x0,xdata,ydata,options);迭代初值迭代初值已知数据点已知数据点选项见无选项见无约束优化约束优化28 6/25/2024 lsqnonlin用以求含参量用以求含参量x x（向量）的向量值函数向量）的向量值函数 f(x)f(x)=(f=(f1 1(x),f(x),f2 2(x),f(x),fn n(x)(x)T T ，使得，使得 最小。最小。其中其中 fi（x）=f（x，xdatai，ydatai）=F(x,xdatai)-ydatai2.lsqnonlin已知数据点：已知数据点：xdataxdata=（xdata1，xdata2，xdataxdatan n）ydataydata=（ydataydata1 1，ydataydata2 2，ydataydatan n）29 6/25/2024 输入格式：输入格式：1）x=lsqnonlin（fun，x0）；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，）；说明：说明：x=lsqnonlinlsqnonlin（fun，x0，options）；）；fun是一个事先建立的是一个事先建立的定义函数定义函数f(x)的的M-文件，文件，自变量为自变量为x迭代初值迭代初值选项见无选项见无约束优化约束优化30 6/25/2024 例例2 用下面一组数据拟合用下面一组数据拟合 中的参数中的参数a，b，k该问题即解的最优化问题：该问题即解的最优化问题：31 6/25/2024 1 1）编写编写M-M-文件文件 curvefun1.mcurvefun1.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）输入命令输入命令tdatatdata=100:100:1000=100:100:1000cdatacdata=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;6.50,6.59;x0=0.2,0.05,0.05;x0=0.2,0.05,0.05;x=x=lsqcurvefitlsqcurvefit(curvefun1,x0,tdata,cdata)(curvefun1,x0,tdata,cdata)f=f=curvefun1(x,tdata)F(x，tdata)=，x=(a，b，k)解法解法1 1.用命令用命令lsqcurvefitlsqcurvefit32 6/25/2024 3 3）运算结果运算结果：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 0.0062 0.0063 0.0063 0.0063 x=0.0063 -0.0034 0.2542 x=0.0063 -0.0034 0.25424）结论）结论：a=0.0063,b=-0.0034,k=0.254233 6/25/2024 1）编写编写M-M-文件文件 curvefun2.mcurvefun2.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)-cdata2）输入命令输入命令:x0=0.2,0.05,0.05;x=lsqnonlin(curvefun2,x0)f=curvefun2(x)函数函数curvefun2的自变量是的自变量是x，cdata和和tdata是已知参数，故应是已知参数，故应将将cdata tdata的值写在的值写在curvefun2.m中中解法解法 2 2 用命令用命令lsqnonlinlsqnonlin x=x=（a a，b b，k k）34 6/25/2024 3 3）运算结果为）运算结果为 f=1.0e-003*f=1.0e-003*（0.2322 -0.1243 -0.2495 -0.2413 0.2322 -0.1243 -0.2495 -0.2413-0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792-0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792）x=0.0063 -0.0034 0.2542x=0.0063 -0.0034 0.2542可以看出，两个命令的计算结果是相同的可以看出，两个命令的计算结果是相同的。4）结论）结论：即拟合得即拟合得a=0.0063 b=-0.0034 k=0.25420.0063 b=-0.0034 k=0.254235 6/25/2024 设有一个年产设有一个年产240吨的食品加工厂吨的食品加工厂,需要统计其生产费用需要统计其生产费用,但由于该厂各项资料不全但由于该厂各项资料不全,无法统计。在这种情况下无法统计。在这种情况下,统统计部门收集了设备、生产能力和该厂大致相同的五个食品计部门收集了设备、生产能力和该厂大致相同的五个食品加工厂的产量与生产费用资料如下表：加工厂的产量与生产费用资料如下表：如何根据已知五个食品加工厂的资料较准确地估计该厂的如何根据已知五个食品加工厂的资料较准确地估计该厂的生产费用呢？生产费用呢？引例引例插值问题插值问题36 6/25/2024 37 6/25/2024 （2）有的函数甚至没有表达式，只是一种表格函数，）有的函数甚至没有表达式，只是一种表格函数，（1）如果函数表达式本身比）如果函数表达式本身比较复复杂，且需要多次，且需要多次实际问题中中经常要涉及到函数常要涉及到函数值的的计算算问题：重复重复计算算时，计算量会很大；算量会很大；方便且表达方便且表达简单的函数来近似代替，的函数来近似代替，这就是就是插插值问题。实际对于于这两种情况，我两种情况，我们都需要都需要寻找一个找一个计算算而我而我们需要的函数需要的函数值可能不在可能不在该表格中。表格中。问题：插：插值与与拟合的区合的区别？38 6/25/2024 一一 维维 插插 值值一、一、插值的定义插值的定义二、插值的方法二、插值的方法三、用三、用Matlab解插值问题解插值问题拉格朗日插值拉格朗日插值牛顿插值牛顿插值反插值反插值三次样条插值三次样条插值分段插值法分段插值法39 6/25/2024 二维插值二维插值一、一、二维插值定义二维插值定义二、二、网格节点插值法网格节点插值法三、三、用用MatlabMatlab解插值问题解插值问题最邻近插值最邻近插值分片线性插值分片线性插值双线性插值双线性插值网格节点数据的插值网格节点数据的插值散点数据的插值散点数据的插值40 6/25/2024 一维插值的定义一维插值的定义已知已知 n+1个节点个节点其中其中互不相同，不妨设互不相同，不妨设求任一插值点求任一插值点处的插值处的插值41 6/25/2024 构造一个构造一个(相对简单的相对简单的)函数函数通过全部节点通过全部节点,即即再用再用计算插值，即计算插值，即42 6/25/2024 代数插值的唯一性代数插值的唯一性是惟一的。是惟一的。定理定理个不同节点，满足插值条件个不同节点，满足插值条件的的n次插值多项式次插值多项式当选择代数多项式作为插值函数时，称为代数多当选择代数多项式作为插值函数时，称为代数多项式插式插值。代数插代数插值43 6/25/2024 一、线性插值一、线性插值（n=1，一次插值），一次插值）求求解解 L1(x)=a1 x+a0已知已知使得使得 L1(xi)=yi.(i=0,1)如果令如果令则称则称 l0(x),l1(x)为为x0,x1上的上的线性插值线性插值基函数基函数。这时。这时根据点斜式得到根据点斜式得到并称其为一次并称其为一次Lagrange插值多项式。插值多项式。f(x)L1(x)=y0l0(x)+y1 l1(x)44 6/25/2024 二、抛物线插值二、抛物线插值（n=2，二次插值），二次插值）求求解解 L2(x)=a2x2+a1 x+a0使得使得 L2(xi)=yi ,i=0,1,2.已知已知45 6/25/2024 抛物插值抛物插值仿照线性函数的构造方法，构造仿照线性函数的构造方法，构造且且其中要求其中要求 均为二次多项式。均为二次多项式。设设即即由由求求46 6/25/2024 故故同理同理47 6/25/2024 例题例题求线性插值函数。求线性插值函数。已知已知 的函数表的函数表 解：解：48 6/25/2024 多项式可使用类似方法。多项式可使用类似方法。分析分析由由 个不同节点个不同节点 构造构造 次次其中其中且且上述多项式称为上述多项式称为拉格朗日插值多项式拉格朗日插值多项式。拉格朗日插值多项式拉格朗日插值多项式49 6/25/2024 插值余项和误差估计插值余项和误差估计且依赖于且依赖于 。其中其中设设f(x)在在a,b有有n+1阶导数，则阶导数，则n次插值多次插值多定理定理项式项式 对任意对任意 ，有插值余项，有插值余项50 6/25/2024 拉格朗日多项式插值的这种振荡现象叫 Runge现象现象 采用拉格朗日多项式插值：选取不同插值节点个数n+1，其中n为插值多项式的次数，当n分别取2,4,6,8,10时，绘出插值结果图形.例例51 6/25/2024 52 6/25/2024 分段线性插值分段线性插值xjxj-1xj+1x0 xnxoy53 6/25/2024 例例用分段线性插值法求插值用分段线性插值法求插值,并观察插值误差并观察插值误差.在在-6,6中平均选取中平均选取41个点作插值个点作插值,结果如图示结果如图示54 6/25/2024 比分段线性插值更光滑。比分段线性插值更光滑。xyxi-1 xiab 在在数数学学上上，光光滑滑程程度度的的定定量量描描述述是是：函函数数(曲曲线线)的的k阶阶导数存在且连续，则称导数存在且连续，则称该曲线具有该曲线具有k阶光滑性阶光滑性。光光滑滑性性的的阶阶次次越越高高，则则越越光光滑滑。是是否否存存在在较较低低次次的的分分段段多多项项式式达达到到较较高高阶阶光光滑滑性性的的方方法法？三三次次样样条条插插值值就就是是一一个很好的例子。个很好的例子。三次样条插值三次样条插值55 6/25/2024 所谓样条函数，从数学上说，就是按照一定光滑性所谓样条函数，从数学上说，就是按照一定光滑性它在每个内节点上具有直到它在每个内节点上具有直到m-1m-1阶连续导数。阶连续导数。要求要求“装配装配”起来的分段多项式。起来的分段多项式。它在每个子区间内是它在每个子区间内是m m次多项式，次多项式，三次样条插值三次样条插值三次样条函数满足：三次样条函数满足：在每个子区间在每个子区间 上是一个三次多项式。上是一个三次多项式。(1)(2)在每个内节点上具有直到二阶的连续导数。在每个内节点上具有直到二阶的连续导数。(3)上可设上可设故在每一个小区间故在每一个小区间求解使用三弯矩法求解使用三弯矩法56 6/25/2024 例例用三次样条插值选取用三次样条插值选取11个基点计算插值个基点计算插值57 6/25/2024 用用MATLABMATLAB作插值计算作插值计算一维插值函数：一维插值函数：yi=interp1(x，y，xi，method)插值方法插值方法被插值点被插值点插值节点插值节点xixi处的插处的插值结果值结果nearest ：最邻近插值最邻近插值linear ：线性插值；线性插值；spline ：三次样条插三次样条插值；值；cubic ：立方插值。立方插值。缺省时：缺省时：分段线性插值。分段线性插值。注意：所有的插值方法都要求注意：所有的插值方法都要求x是单调的，并且是单调的，并且xi不能不能够超过够超过x的范围。的范围。58 6/25/2024 例：在例：在1-121-12的的1111小时内，每隔小时内，每隔1 1小时测量一次温小时测量一次温度，测得的温度依次为：度，测得的温度依次为：5 5，8 8，9 9，1515，2525，2929，3131，3030，2222，2525，2727，2424。试估计每隔。试估计每隔1/101/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)59 6/25/2024 60 6/25/2024 xy机翼下轮廓线例例 已知飞机下轮廓线上数据如下，求已知飞机下轮廓线上数据如下，求x每改变每改变0.1时的时的y值。值。61 6/25/2024 62 6/25/2024 二维插值的定义二维插值的定义 xyO O第一种（网格节点）：第一种（网格节点）：63 6/25/2024 已知已知 m n个节点个节点 其中其中互不相同，不妨设互不相同，不妨设 构造一个二元函数构造一个二元函数通过全部已知节点通过全部已知节点,即即再用再用计算插值，即计算插值，即64 6/25/2024 第二种（散乱节点）：第二种（散乱节点）：yx0 065 6/25/2024 已知已知n个节点个节点其中其中互不相同，互不相同，构造一个二元函数构造一个二元函数通过全部已知节点通过全部已知节点,即即再用再用计算插值，即计算插值，即66 6/25/2024 注意：注意：最邻近插值一般不连续。具有连续性的最简单的插值是分片线性插值。最邻近插值最邻近插值x y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O O 二维或高维情形的最邻近插值，与被插值点最邻近的节点的函数值即为所求。67 6/25/2024 将四个插值点（矩形的四个顶点）处的函数值依次简记为：分片线性插值分片线性插值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)=f468 6/25/2024 插值函数为：第二片(上三角形区域)：(x,y)满足插值函数为：注意注意：(x,y)当然应该是在插值节点所形成的矩形区域内。显然，分片线性插值函数是连续的；分两片的函数表达式如下：第一片(下三角形区域)：(x,y)满足69 6/25/2024 双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下：其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。双线性插值双线性插值x y(x1,y1)(x1,y2)(x2,y1)(x2,y2)O O70 6/25/2024 要求要求x0,y0 x0,y0单调；单调；x x，y y可取可取为矩阵，或为矩阵，或x x取取行向量，行向量，y y取为列向量，取为列向量，x,yx,y的值分别不能超出的值分别不能超出x0,y0 x0,y0的范围。的范围。z=interp2(x0,y0,z0,x,y,method)被插值点插值方法用用MATLAB作网格节点数据的插值作网格节点数据的插值插值节点被插值点的函数值nearest nearest 最邻近插值最邻近插值linear linear 双线性插值双线性插值cubic cubic 双三次插值双三次插值缺省时缺省时,双线性插值双线性插值71 6/25/2024 例：测得平板表面例：测得平板表面3*53*5网格点处的温度分别为：网格点处的温度分别为：82 81 80 82 84 82 81 80 82 84 79 63 61 65 81 79 63 61 65 81 84 84 8484 82 85 86 82 85 86 试作出平板表面的温度分布曲面试作出平板表面的温度分布曲面z=f(x,y)z=f(x,y)的图形。的图形。输入以下命令：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.先在三维坐标画出原始数据，画出粗糙的温度分布曲图.72 6/25/2024 73 6/25/2024 再输入以下命令: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方向上每隔0.2个单位的地方进行插值.74 6/25/2024 75 6/25/2024 通过此例对最近邻点插值、双线性插值方法和双三次插值方法的插值效果进行比较。76 6/25/2024 原始数据图原始数据图77 6/25/2024 最邻近插值最邻近插值78 6/25/2024 双线性插值79 6/25/2024 双三次插值80 6/25/2024 插值函数插值函数griddata格式为格式为:cz=griddata（x，y，z，cx，cy，method）用用MATLABMATLAB作散点数据的插值计算作散点数据的插值计算 要求要求cxcx取行向量，取行向量，cycy取为列向量取为列向量。被插值点插值方法插值节点被插值点的函数值nearest nearest 最邻近插值最邻近插值linear linear 双线性插值双线性插值cubic cubic 双三次插值双三次插值v4-Matlab提供的插值方法提供的插值方法缺省时缺省时,双线性插值双线性插值81 6/25/2024 一室模型一室模型：将整个机体看作一个房室，称：将整个机体看作一个房室，称中心室中心室，室内血药浓度是均匀的。快速静脉注射后，浓度立室内血药浓度是均匀的。快速静脉注射后，浓度立即上升；然后迅速下降。当浓度太低时，达不到预即上升；然后迅速下降。当浓度太低时，达不到预期的治疗效果；当浓度太高，又可能导致药物中毒期的治疗效果；当浓度太高，又可能导致药物中毒或副作用太强。临床上，每种药物有一个最小有效或副作用太强。临床上，每种药物有一个最小有效浓度浓度c1和一个最大有效浓度和一个最大有效浓度c2。设计给药方案时，要设计给药方案时，要使血药浓度使血药浓度 保持在保持在c1c2之间。本题设之间。本题设c1=10，c2=25(ug/ml).拟拟 合合 问问 题题 实实 例例给药方案给药方案 一种新药用于临床之前，必须设计给药方案一种新药用于临床之前，必须设计给药方案.药物进入机体后血液输送到全身，在这个过程药物进入机体后血液输送到全身，在这个过程中不断地被吸收、分布、代谢，最终排出体外，药中不断地被吸收、分布、代谢，最终排出体外，药物在血液中的浓度，即单位体积血液中的药物含量，物在血液中的浓度，即单位体积血液中的药物含量，称为称为血药浓度。血药浓度。82 6/25/2024 在实验方面在实验方面,对某人用快速静脉注射方式一次对某人用快速静脉注射方式一次注入该药物注入该药物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 要设计给药方案要设计给药方案,必须知道给药后血药浓度随必须知道给药后血药浓度随时间变化的规律。从实验和理论两方面着手：时间变化的规律。从实验和理论两方面着手：83 6/25/2024 给药方案给药方案 1.在快速静脉注射的给药方式下，研究血药浓在快速静脉注射的给药方式下，研究血药浓度（单位体积血液中的药物含量）的变化规律。度（单位体积血液中的药物含量）的变化规律。tc2cc10问问题题2.给定药物的最小有效浓度和最大治疗浓度，设给定药物的最小有效浓度和最大治疗浓度，设计给药方案：每次注射剂量多大；间隔时间多长。计给药方案：每次注射剂量多大；间隔时间多长。分分析析 理论：用一室模型研理论：用一室模型研究血药浓度变化规律究血药浓度变化规律 实验：对血药浓实验：对血药浓度数据作拟合，符度数据作拟合，符合负指数变化规律合负指数变化规律84 6/25/2024 3.3.血液容积血液容积v,t=0v,t=0注射剂量注射剂量d,d,血药浓度立即为血药浓度立即为d/v.d/v.2.2.药物排除速率与血药浓度成正比，比例系数药物排除速率与血药浓度成正比，比例系数 k(0)k(0)模型假设模型假设1.1.机体看作一个房室，室内血药浓度均匀机体看作一个房室，室内血药浓度均匀一室模型一室模型模型建立模型建立 在此，在此，d=300mg，t及及c（t）在某些点处的值见前表，在某些点处的值见前表，需经拟合求出参数需经拟合求出参数k、v85 6/25/2024 用线性最小二乘拟合用线性最小二乘拟合c(t)计算结果：计算结果：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)程序：程序：86 6/25/2024 给药方案给药方案 设计设计cc2c10t 设每次注射剂量D,间隔时间 血药浓度c(t)应c1 c(t)c2 初次剂量D0 应加大给药方案记为：给药方案记为：2、1、计算结果：计算结果：给药方案：给药方案：c1=10,c2=25k=0.2347v=15.0287 6/25/2024 故可制定给药方案：故可制定给药方案：即即:首次注射首次注射375mg，其余每次注射其余每次注射225mg，注射的间隔时间为注射的间隔时间为4小时。小时。