教学项目 数值法模型

数学建模与实验严可颂主讲教学项目四数值分析法模型4.1插值法建模

拉格朗日插值分段线性插值三次样条插值一、插值的定义二、插值的方法三、用Matlab解插值问题已知n+1个节点其中互不相同，不妨设求任一插值点处的插值节点可视为由产生，表达式复杂，或无封闭形式，或未知。4.1.1插值法的定义

构造一个(相对简单的)函数通过全部节点,即再用计算插值，即4.1.2插值的方法

已知函数f(x)在n+1个点x0,x1,…,xn处的函数值为y0,y1,…,yn

。求一n次多项式函数Pn(x)，使其满足：拉格朗日(Lagrange)插值称为拉格朗日插值基函数。解决此问题的拉格朗日插值多项式公式如下其中Li(x)为n次多项式：特殊情况两点一次(线性)插值多项式:三点二次(抛物)插值多项式:采用拉格朗日多项式插值：选取不同插值节点个数n+1，其中n为插值多项式的次数，当n取10时，绘出插值结果图形.例分段线性插值计算量与n无关;n越大，误差越小.xjxj-1xj+1x0xnxoy比分段线线性插值值更光滑滑。xyxi-1xiab在数学上上，光滑滑程度的的定量描描述是：：函数(曲线)的k阶导数存存在且连连续，则则称该曲线具具有k阶光滑性性。光滑性的的阶次越越高，则则越光滑滑。是否否存在较较低次的的分段多多项式达达到较高高阶光滑滑性的方方法？三三次样条条插值就就是一个个很好的的例子。。三次样条条插值三次样条条插值g(x)为被插值函函数。用MATLAB作插值计计算一维插值值函数：：yi=interp1(x，y，xi，'method')插值方法被插值点插值节点xi处的插值结果‘nearest’：最邻近近插值‘linear’’：线性插值值；‘spline’’：三次次样条插插值；‘cubic’：立方方插值。。缺省时：：分分段线性性插值。。注意：所所有的插插值方法法都要求求x是单调的的，并且且xi不能够超超过x的范围。。例：在1-12的11小时内，，每隔1小时测量量一次温温度，测测得的温温度依次次为：5，8，9，15，25，29，31，30，22，25，27，24。试估计计每隔1/10小时的温温度值。。hours=1:12;h=1:0.1:12;t=interp1(hours,temps,h,'spline');(直接输出出数据将将是很多多的)plot(hours,temps,'+',h,t,hours,temps,'r:')%作图xlabel('Hour'),ylabel('DegreesCelsius’)xy实验习题题已知飞机机下轮廓廓线上数数据如下下，求x每改变0.1时的y值。机翼下轮廓线4.2拟合法简介

2.拟合的基基本原理理1.拟合问题题引例拟合问问题题引例例1温度t(0C)20.532.751.073.095.7电阻R()7658268739421032已知热敏电阻数据：求600C时的电阻阻R。设R=at+ba,b为待定系系数拟合问问题题引例例2

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01已知一室模型快速静脉注射下的血药浓度数据(t=0注射300mg)求血药浓浓度随时时间的变变化规律律c(t).作半对数数坐标系系(semilogy)下的图形形曲线拟拟合合问题题的的提法法已知一组组（二维维）数据据，即平平面上n个点（xi,yi)i=1,…n,寻求一个个函数（（曲线））y=f(x),使f(x)在某种准准则下与与所有数数据点最最为接近近，即曲曲线拟合合得最好好。+++++++++xyy=f(x)(xi,yi)ii为点（xi,yi)与曲线y=f(x)的距离拟合与插值的的关系函数插值与曲曲线拟合都是是要根据一组组数据构造一一个函数作为为近似，由于于近似的要求求不同，二者者的数学方法法上是完全不不同的。实例：下面数据是某某次实验所得得，希望得到到X和f之间的关系？？问题：给定一批数据据点，需确定定满足特定要要求的曲线或或曲面解决方案：若不要求曲线线（面）通过过所有数据点点，而是要求求它反映对象象整体的变化化趋势，这就就是数据拟合，又称曲线拟拟合或曲面拟拟合。若要求所求曲曲线（面）通通过所给所有有数据点，就就是插值问题；最临近插值、、线性插值、、样条插值与与曲线拟合结结果：曲线拟合问题题最常用的解解法——线性最小二乘乘法的基本思思路第一步:先选定一组函函数r1(x),r2(x),……rm(x),m<n,令f(x)=a1r1(x)+a2r2(x)+……+amrm(x)（1）其中a1,a2,…am为待定系数。。第二步:确定a1,a2,…am的准则（最小小二乘准则））：使n个点（xi,yi)与曲线y=f(x)的距离i的平方和最小小。记

问题归结为，，求a1,a2,…am使J(a1,a2,…am)最小。线性最小二乘乘法的求解：：预备知识超定方程组：方程个数大大于未知量个个数的方程组组即Ra=y其中超定方程一般般是不存在解解的矛盾方程程组。如果有向量a使得达到最小，则称a为上述超定方程的最小二乘解。线性最小二乘乘法的求解定理：当RTR可逆时，超定定方程组（3）存在最小二二乘解，且即即为方程组RTRa=RTy的解：a=(RTR)-1RTy所以，曲线拟拟合的最小二二乘法要解决决的问题，实实际上就是求求以下超定方方程组的最小小二乘解的问问题。其中Ra=y（3）线性最小二乘乘拟合f(x)=a1r1(x)+……+amrm(x)中函数{r1(x),……rm(x)}的选取1.通过机理分析析建立数学模模型来确定f(x)；++++++++++++++++++++++++++++++f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx2.将数据(xi,yi)i=1,…n作图，通过直直观判断确定定f(x)：用MATLAB解拟合问题1、线性最小二二乘拟合2、非线性最小小二乘拟合用MATLAB作线性最小二二乘拟合1.作多项式f(x)=a1xm+…+amx+am+1拟合,可利用已有程程序:a=polyfit(x,y,m)2.对超定方程组可得最小二乘意义下的解。，用3.多项式在x处的值y可用以下命令令计算：y=polyval（a，x）输入同长度的数组X，Y拟合多项式次数即要求出二次多项式:中的使得:例对下面一组数据作二次多项式拟合1）输入以下命命令：x=0:0.1:1;y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2];R=[(x.^2)'x'ones(11,1)]；A=R\y'解法1．用解超定方程程的方法2）计算结果：Ａ1）输入以下命命令：x=0:0.1:1;y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2];A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,'k+',x,z,'r')%作出数据点和和拟合曲线的的图形2）计算结果：：Ａ解法2．用多项式拟合合的命令1.lsqcurvefit已知数据点：xdata=（xdata1，xdata2，…，xdatan），ydata=（ydata1，ydata2，…，ydatan）用MATLAB作非线性最小小二乘拟合Matlab的提供了两个个求非线性最最小二乘拟合合的函数：lsqcurvefit和lsqnonlin。两个命令都都要先建立M-文件fun.m，在其中定义义函数f(x)，但两者定义义f(x)的方式是不同同的,可参考例题.

lsqcurvefit用以求含参量x（向量）的向量值函数F(x,xdata)=（F（x，xdata1），…，F（x，xdatan））T中的参变量x(向量),使得输入格式为:(1)x=lsqcurvefit(‘fun’’,x0,xdata,ydata);(2)x=lsqcurvefit(‘fun’’,x0,xdata,ydata,options);(3)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options,’grad’);(4)[x,options]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(5)[x,options,funval]=lsqcurvefit(‘‘fun’,x0,xdata,ydata,…);(6)[x,options,funval,Jacob]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);fun是一个事先建立的定义函数F(x,xdata)

的M-文件,自变量为x和xdata说明：x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);迭代初值已知数据点选项见无约束优化

lsqnonlin用以求含参量x（向量）的向量值函数

f(x)=(f1(x),f2(x),…,fn(x))T

中的参量x，使得

最小。其中fi（x）=f（x，xdatai，ydatai）

=F(x,xdatai)-ydatai

2.lsqnonlin已知数据点：：xdata=（xdata1，xdata2，…，xdatan）ydata=（ydata1，ydata2，…，ydatan）输入格式为：：1）x=lsqnonlin（‘fun’，x0）；2）x=lsqnonlin（‘fun’，x0，options）；3）x=lsqnonlin（‘fun’，x0，options，‘grad’）；4）[x，options]=lsqnonlin（‘fun’，x0，…）；5）[x，options，funval]=lsqnonlin（‘fun’，x0，…）；说明：x=lsqnonlin（‘fun’，x0，options）；fun是一个事先建立的定义函数f(x)的M-文件，自变量为x迭代初值选项见无约束优化

例2用下面一组数据拟合中的参数a，b，k该问题即解最最优化问题：：1）编写M-文件curvefun1.mfunctionf=curvefun1(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata)%其中x(1)=a;x(2)=b；x(3)=k;2）输入命令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,0.05];x=lsqcurvefit('curvefun1',x0,tdata,cdata)f=curvefun1(x,tdata)

F(x，tdata)=，x=(a，b，k)解法1.用命令令lsqcurvefit3）运算算结果果为：4）结论论：a=0.0063,b=-0.0034,k=0.2542

解法2

用命令lsqnonlin

f(x)=F(x,tdata,ctada)=x=（a，b，k）1）编写M-文件curvefun2.mfunctionf=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)3）运算算结果果为可以看看出,两个命命令的的计算算结果果是相相同的的.4）结论论：即拟合合得a=0.0063b=-0.0034k=0.2542MATLAB解应用用问题题实例例1、电阻阻问题题2、给药药方案案问题题*3、水塔塔流量量估计计问题题电阻问问题温度t(0C)20.532.751.073.095.7电阻R()7658268739421032例.由数据拟合R=a1t+a2方法1.用命令polyfit(x,y,m)得到a1=3.3940,a2=702.4918方法2.直接用结果相相同。。一室模模型：将整整个机机体看看作一一个房房室，，称中心室室，室内内血药药浓度度是均均匀的的。快快速静静脉注注射后后，浓浓度立立即上上升；；然后后迅速速下降降。当当浓度度太低低时，，达不不到预预期的的治疗疗效果果；当当浓度度太高高，又又可能能导致致药物物中毒毒或副副作用用太强强。临临床上上，每每种药药物有有一个个最小小有效效浓度度c1和一个个最大大有效效浓度度c2。设计计给药药方案案时，，要使使血药药浓度度保保持在在c1~c2之间。。本题题设c1=10，c2=25(ug/ml).拟合问题实例2给药方案——一种新新药用用于临临床之之前，，必须须设计计给药药方案案.药物进进入机机体后后血液液输送送到全全身，，在这这个过过程中中不断断地被被吸收收、分分布、、代谢谢，最最终排排出体体外，，药物物在血血液中中的浓浓度，，即单单位体体积血血液中中的药药物含含量，，称为为血药浓浓度。。

在实验方面,对某人用快速静脉注射方式一次注入该药物300mg后,在一定时刻t(小时)采集血药,测得血药浓度c(ug/ml)如下表:

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01要设计计给药药方案案,必须知知道给给药后后血药药浓度度随时时间变变化的的规律律。从从实验验和理理论两两方面面着手手：给药方方案1.在快速速静脉脉注射射的给给药方方式下下，研研究血血药浓浓度（（单位位体积积血液液中的的药物物含量量）的的变化化规律律。tc2cc10问题题2.给定定药药物物的的最最小小有有效效浓浓度度和和最最大大治治疗疗浓浓度度，，设设计计给给药药方方案案：：每每次次注注射射剂剂量量多多大大；；间间隔隔时时间间多多长长。。分析析理论论：：用用一一室室模模型型研研究究血血药药浓浓度度变变化化规规律律实验验：：对对血血药药浓浓度度数数据据作作拟拟合合，，符符合合负负指指数数变变化化规规律律3.血液液容容积积v,t=0注射射剂剂量量d,血药药浓浓度度立立即即为为d/v.2.药物物排排除除速速率率与与血血药药浓浓度度成成正正比比，，比比例例系系数数k(>0)模型型假假设设1.机体体看看作作一一个个房房室室，，室室内内血血药药浓浓度度均均匀匀———一室室模模型型模型型建建立立在此此，，d=300mg，t及c（t）在在某某些些点点处处的的值值见见前前表表，，需需经经拟拟合合求求出出参参数数k、v用线线性性最最小小二二乘乘拟拟合合c(t)计算结果：d=300;t=[0.250.511.523468];c=[19.2118.1515.3614.1012.899.327.455.243.01];y=log(c);a=polyfit(t,y,1)k=-a(1)v=d/exp(a(2))程序：给药药方方案案设设计计cc2c10t设每每次次注注射射剂剂量量D,间隔隔时时间间血药药浓浓度度c(t)应c1c(t)c2初次次剂剂量量D0应加加大大给药方案记为：2、1、计算结果：给药方案：c1=10,c2=25k=0.2347v=15.02故可可制制定定给给药药方方案案：：即:首次次注注射射375mg，其余余每每次次注注射射225mg，注射射的的间间隔隔时时间间为为4小时时。。估计计水水塔塔的的流流量量2、解解题题思思路路3、算算法法设设计计与与编编程程1、问问题题某居居民民区区有有一一供供居居民民用用水水的的园园柱柱形形水水塔塔，，一一般般可可以以通通过过测测量量其其水水位位来来估估计计水水的的流流量量，，但但面面临临的的困困难难是是，，当当水水塔塔水水位位下下降降到到设设定定的的最最低低水水位位时时，，水水泵泵自自动动启启动动向向水水塔塔供供水水，，到到设设定定的的最最高高水水位位时时停停止止供供水水，，这这段段时时间间无无法法测测量量水水塔塔的的水水位位和和水水泵泵的的供供水水量量．．通通常常水水泵泵每每天天供供水水一一两两次次，，每每次次约约两两小小时时.水塔塔是是一一个个高高12.2米，，直直径径17.4米的的正正园园柱柱．．按按照照设设计计，，水水塔塔水水位位降降至至约约8.2米时时，，水水泵泵自自动动启启动动，，水水位位升升到到约约10.8米时时水水泵泵停停止止工工作作．．表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、2时段的水水位，并并导出流流量2、拟合供水水时段的的流量3、估计一天天总用水水量4、流量及及总用水水量的检检验1、拟合第1时段的水水位，并并导出流流量设t，h为已输入入的时刻刻和水位位测量记记录（水水泵启动动的4个时刻不不输入）），第1时段各时刻的的流量可可如下得得：1）c1=polyfit（t（1：10），h（1：10），3）；%用3次多项式式拟合第第1时段水位位，c1输出3次多项式式的系数数2）a1=polyder（c1）；%a1输出多项项式（系系数为c1）导数的的系数3）tp1=0：0.1：9；x1=-polyval（a1，tp1）；%x1输出多项项式（系系数为a1）在tp1点的函数数值（取取负后边边为正值值），即即tp1时刻的流流量4）流量函数数为：2、拟合第2时段的水水位，并并导出流流量设t，h为已输入入的时刻刻和水位位测量记记录（水水泵启动动的4个时刻不不输入）），第2时段各时刻的的流量可可如下得得：1）c2=polyfit(t(10.9:21),h(10.9:21),3)；%用3次多项式式拟合第第2时段水位位，c2输出3次多项式式的系数数2）a2=polyder(c2)；%a2输出多项项式（系系数为c2）导数的的系数3）tp2=10.9:0.1:21;x2=-polyval(a2,tp2);%x2输出多项项式（系系数为a2）在tp2点的函数数值（取取负后边边为正值值），即即tp2时刻的流流量4）流量函数数为：3、拟合供水水时段的的流量在第1供水时段段（t=9~11）之前（（即第1时段）和和之后（（即第2时段）各各取几点点，其流流量已经经得到，，用它们们拟合第第1供水时段段的流量量．为使使流量函函数在t=9和t=11连续，我我们简单单地只取取4个点，拟拟合3次多项式式（即曲曲线必过过这4个点），，实现如如下：xx1=-polyval（a1，[89]）；%取第1时段在t=8,9的流量xx2=-polyval（a2，[1112]）；%取第2时段在t=11,12的流量xx12=[xx1xx2]；c12=polyfit（[891112]，xx12，3）；%拟合3次多项式式tp12=9：0.1：11；x12=polyval（c12，tp12）；%x12输出第1供水时段段各时刻的的流量拟合的流流量函数数为：在第2供水时段段之前取取t=20，20.8两点的流流水量，，在该时时刻之后后（第3时段）仅仅有3个水位记记录，我我们用差差分得到