数据拟合.docx

1、数据拟合数据拟合根据一组二维数据，即平面上的若干点，要求确定一个一元函数y = f(x)，即曲线，使这些点与曲线总体来说尽量接近。这就是数据拟合成曲线的思想，简称为曲线拟合(fitting a curve)。曲线拟合其目的是根据实验获得的数据去建立因变量与自变量之间有效的经验函数关系，为进一步的深入研究提供线索。本章的目的，掌握一些曲线拟合的基本方法，弄清楚曲线拟合与插值方法之间的区别,学会使用MATLAB软件进行曲线拟合。5.1 引例拟合问题引例一 电阻问题 已知热敏电阻电阻值与温度的数据：温度t(0C)20.532.751.073.095.7电阻R()7658268739421032求温度

2、为63度时的电阻值。拟合问题引例二 给药问题 一种新药用于临床之前，必须设计给药方案。 药物进入机体后血液输送到全身，在这个过程中不断地被吸收、分布、代谢，最终排出体外，药物在血液中的浓度，即单位体积血液中的药物含量，称为血药浓度。 一室模型：将整个机体看作一个房室，称中心室， 室内血药浓度是均匀的。快速静脉注射后，浓度立即上升；然后迅速下降。 当浓度太低时，达不到预期的治疗效果；当浓度太高，又可能 导致药物中毒或副作用太强。 临床上，每种药物有一个最小有效浓度c1和一个最大有效浓度c2。设计给药方案时，要使血药浓度 保持在c1c2之间。 本题设c1=10，c2=25(ug/ml). 要设计给

3、药方案,必须知道给药后血药浓度随时间变化的规律。从实验和理论两方面着手：在实验方面, t=0时对某人用快速静脉注射方式一次注入该药物300mg后,在一定时刻t(小时)采集血药,测得血药浓度c(ug/ml)如下表:t (h)0.250.511.523468c (g/ml)19.2118.1515.3614.1012.89 9.327.455.24 3.01问题：1. 在快速静脉注射的给药方式下，研究血药浓度（单位体积血液中的药物含量）的变化规律。2. 给定药物的最小有效浓度和最大治疗浓度，设计给药方案：每次注射剂量多大；间隔时间多长。5.2 最小二乘法给定平面上的点(xi, yi)，（i = 1

4、,2,n），进行曲线拟合有多种方法，其中最小二乘法是解决曲线拟合最常用的方法。最小二乘法的原理是：求 达到最小。如图1所示，其中i为点(xi,yi)与曲线y=f(x)的距离。曲线拟合的实际含义是寻求一个函数y=f(x)，使f(x)在某种准则下与所有数据点最为接近，即曲线拟合得最好。最小二乘准则就是使所有散点到曲线的距离平方和最小。拟合时选用一定的拟合函数f(x) 形式，设拟合函数可由一些简单的“基函数”（例如幂函数，三角函数等等） 来线性表示：图1 曲线拟合示意图现在要确定系数c0，c1，cm，使达到极小。为此，将f(x)的表达式代入中，就成为c0，c1，cm的函数，求的极小，就可令对 ci的

5、偏导数等于零，于是得到m+1个方程组，从中求解出ci。通常取基函数为1，x，x2，x3，xm，这时拟合函数f(x)为多项式函数。当m=1时，f(x) = a + bx，称为一元线性拟合函数，它是曲线拟合最简单的形式。除此之外，常用的一元曲线拟合函数还有双曲线f(x) = a + b/x，指数曲线f(x) = a e bx 等，对于这些曲线，拟合前须作变量代换，转化为线性函数。已知一组数据，用什么样的曲线拟合最好呢？可以根据散点图进行直观判断，在此基础上，选择几种曲线分别作拟合，然后比较，观察哪条曲线的最小二乘指标最小。5.3 曲线拟合的MATLAB实现MATLAB软件提供了基本的曲线拟合函数的

6、命令：1多项式函数拟合： a = polyfit（xdata，ydata，n）其中n表示多项式的最高阶数，xdata，ydata 为要拟合的数据，它是用数组的方式输入。输出参数a为拟合多项式y = a1xn + + anx + an+1的系数a = a1, , an, an+1。多项式在x处的值y可用下面程序计算。y = polyval (a, x)2一般的曲线拟合： p = curvefit（Funp0，xdata，ydata）其中Fun表示函数Fun (p, xdata)的M-文件，p0表示函数的初值。curvefit命令的求解问题形式是：minp sum (Fun (p, xdata)-

7、ydata).2若要求解点x处的函数值可用程序f = Fun(p, x) 计算。例如已知函数形式，并且已知数据点(xi, yi), i = 1,2, n，要确定四个未知参数a, b, c, d。使用curvefit命令，数据输入xdata = x1,x2, , xn; ydata = y1,y2, , yn;初值输入p0 = a0,b0,c0,d0; 并且建立函数的M-文件（Fun.m）。若定义p1 = a, p2 = b, p3 = c, p4 = d , 则输出p = p1, p2, p3, p4。3引例1求解 根据热敏电阻电阻值与温度的数据：温度t(0C)20.532.751.073.0

8、95.7电阻R()7658268739421032首先作出散点图t=20.5 32.7 51.0 73.0 95.7;r=765 826 873 942 1032;plot(t,r) 从散点图看出电阻值和温度之间的关系近似于线性。因此设R=a1 t +a2 a=polyfit(t,r,1) a = 3.3987 702.0968 表示R=3.3987 t +702.096863度时的电阻值为r=polyval(a,63) r = 916.2174 4引例2求解：（1）模型假设1. 机体看作一个房室，室内血药浓度均匀一室模型2. 药物排除速率与血药浓度c成正比，比例系数 k(0)3. 血液容积

9、v, t=0时注射剂量 d, 血药浓度即为 d/v.（2）模型建立由假设2得： 由假设3得： 在此，d=300mg，t及c（t）在某些点处的值见前表，需经拟合求出参数k、vd=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)计算结果：a = -0.2347 2.9943k = 0.2347v = 15.0219画出血药浓度随时间的变化规律图t1=0:0.1:8;ct=(d/v)*exp(-k

10、*t1);plot(t,c,o,t1,ct,g-) 设每次注射剂量D, 间隔时间 血药浓度c(t) 应c1 c(t) c2 初次剂量D0 应加大给药方案记为： 1. 2. 其中c1=10，c2=25，k=0.2347，v=15.02计算结果：给药方案：即: 首次注射 375 mg， 其余每次注射 225 mg， 注射的间隔时间为 4 小时。 5.拟合与插值的比较问题：给定一批离散的数据点，需确定满足特定要求的曲线或曲面, 从而获取整体的规律。即通过窥几斑来达到知全豹。 解决方案： 若要求所求曲线（面）通过所给所有数据点，就是插值问题； 若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的

11、变化趋势，这就是数据拟合，又称曲线拟合或曲面拟合。 从几何意义上看，拟合是给定了空间中的一些点，找到一个已知形式的连续曲面来最大限度地逼近这些点；而插值是找到一个(或几个分片光滑的)连续曲面来穿过这些点。 附：书上引例：浓度变化规律的计算x=1:16;y=4,6.4,8,8.4,9.28,9.5,9.7,9.86,10,10.2,10.32,10.42,10.5,10.55,10.58,10.6; plot(x,y,o) p=polyfit(x,y,2) 计算结果： p = -0.0445 1.0711 4.3252 %二次多项式的系数从而得到某化合物的浓度y与时间x的拟合函数：对函数的精度如

12、何检测呢？仍然以图形来检测，将散点与拟合曲线画在一个画面上。 xi=linspace(0,16,160); yi=polyval(p,xi); plot(x,y,o,xi,yi) （这两段程序在此无法运行，需在MATLAB的命令窗口中运行）参见图2。由此看出上述曲线拟合是比较吻合的。图2 浓度y的拟合曲线与实测数据(o)的比较下面我们来看看曲线拟合与曲线插值有什么区别：MATLAB程序如下：t=interp1(x,y,xi,linear);plot(x,y,+,xi,t, r:)hold onplot(xi,yi) 该程序运行后，得到图3所示的结果。其中，标有+的是已知数据点，连接数据点的虚线

13、是线性插值函数曲线，光滑的函数曲线是最佳拟合曲线。图3 拟合曲线与线性插值曲线的比较在MATLAB的NAG Foundation Toolbox中也有一些曲面拟合函数，如e02daf，e02cf，e02def可分别求出矩形网格点数据、散点数据的最小平方误差双三次样条曲面拟合，e02def等可求出曲面拟合的函数值。5.4 范例：薄膜渗透率的测定某种医用薄膜有允许一种物质的分子穿透它，从高浓度的溶液向低浓度的溶液扩散的功能，在试制时需测定薄膜被这种分子穿透的能力。测定方法如下：用面积为S的薄膜将容器分成体积分别为的两部分，在两部分中分别注满该物质的两种不同浓度的溶液。此时该物质分子就会从高浓度溶液

14、穿过薄膜向低浓度溶液中扩散。通过单位面积薄膜分子扩散的速度与膜两侧溶液的浓度差成正比，比例系数K表征了薄膜被该物质分子穿透的能力，称为渗透率。定时测量容器中薄膜某一侧的溶液浓度值，以此确定K的数值。541 数学模型的建立这是一个综合性的应用实例。主要涉及微分方程和数据拟合参数（即参数辨识）的数学知识。1假设：1） 薄膜两侧的溶液始终是均匀的，即在任何时刻膜两侧的每一处溶液的浓度都是相同的。2） 当两侧浓度不一致时，物质的分子穿透薄膜总是从高浓度溶液向低浓度溶液扩散。3） 通过单位面积膜分子扩散的速度与膜两侧溶液的浓度差成正比。4） 薄膜是双向同性的即物质从膜的任何一侧向另一侧渗透的性能是相同的

15、。图4 圆柱体容器被薄膜截面S阻隔不同的假设应该具有不同形式的数学模型。2符号说明： 1）表示t时刻膜两侧溶液的浓度;2）表示初始时刻两侧溶液的浓度 (单位:毫克/立方厘米);3）K 表示渗透率;4）表示由薄膜阻隔的容器两侧的体积;3分析：考察时段t，t+t薄膜两侧容器中该物质质量的变化。以容器A侧为例，在该时段物质质量的增加量： 另一方面由渗透率的定义我们知道，从B侧渗透至A侧的该物质的质量为：。由质量守恒定律，两者应该相等，于是有=。两边除以t，令t0并整理得（5.1）且注意到整个容器的溶液中含有该物质的质量应该不变，即有下式成立： 代入（5.1）得： 再利用初始条件 解出： 至此，问题归

16、结为利用来辨识参数K和，对应的数学模型变为求函数： 令 问题转化为求函数 的最小值点（K，a，b）。5.4.2 求解参数例如，设立方厘米，S =10平方厘米，对容器的B部分溶液浓度的测试结果如表1。表1tj (s)1002003004005006007008009001000Cj(10-5)4.544.995.355.655.906.106.266.396.506.59其中Cj的单位：毫克/立方厘米。此时极小化的函数为: 用MATLAB软件进行计算。1） 编写M-文件 curvefun.m function f = curvefun(x,tdata) f = x(1)+x(2)*exp(-0.

17、02*x(3)*tdata);其中 x(1) = a;x(2) = b;x(3) = k;2）编写程序(test1.m) tdata = linspace(100,1000,10); cdata = 1e-05.*454 499 535 565 590 610 626 639 650 659; x0 = 0.2,0.05,0.05; x = curvefit(curvefun,x0,tdata,cdata) f= curvefun(x,tdata)e=f-cdataplot(tdata,cdata,o,tdata,f,r-) 3) 输出结果:x = 0.0070 -0.0030 0.1012即

18、表示k = 0.1012, a = 0.007, b = -0.003f = Columns 1 through 8 0.0045 0.0050 0.0054 0.0057 0.0059 0.0061 0.0063 0.0064 Columns 9 through 10 1.65 0.0066误差值为：e = 1.0e-005 * Columns 1 through 8 -0.0339 -0.2146 0.3900 0.2857 -0.2981 -0.3570 -0.0707 0.2305 Columns 9 through 10 0.0938 -0.0339现在curvefit函数已经被ls

19、qcurvefit函数取代tdata = linspace(100,1000,10); cdata = 1e-05.*454 499 535 565 590 610 626 639 650 659; x0 = 0.2,0.05,0.05; x = lsqcurvefit(curvefun,x0,tdata,cdata) f= curvefun(x,tdata)e=f-cdataplot(tdata,cdata,o,tdata,f,r-) 运行结果x = 0.0063 -0.0034 0.2542f = Columns 1 through 8 0.0043 0.0051 0.0056 0.005

20、9 0.0061 0.0062 0.0062 0.0063 Columns 9 through 10 1.63 0.0063e = 1.0e-003 * Columns 1 through 8 -0.2322 0.1243 0.2495 0.2413 0.1668 0.0724 -0.0241 -0.1159 Columns 9 through 10 -0.2030 -0.2792 (该程序在MATLAB命令窗口运行)进一步求得：（毫克/立方厘米），（毫克/立方厘米）在MATLAB中,还有一种求解非线性数据拟合的语句。标准问题: min f(x) = f1(x)2 + f2(x)2 + +fm

21、(x)2 + L 其中 L为常数。基本格式:x = leastsq (fun,x0) 其中fun为fi(x)(i = 1,2,m)函数的M-文件。5.5 数据插值与拟合实验一、实验目的及意义1 了解插值、最小二乘拟合的基本原理2 掌握用MATLAB计算一维插值和两种二维插值的方法；3 掌握用MATLAB作最小二乘多项式拟合和曲线拟合的方法。二、实验内容1针对实际问题，试建立数学模型。用MATLAB计算一维插值和两种二维插值的方法求解；1用MATLAB中的函数作一元函数的多项式拟合与曲线拟合，作出误差图；2用MATLAB中的函数作二元函数的最小二乘拟合，作出误差图；3针对预测和确定参数的实际问题

22、，建立数学模型，并求解。三、实验步骤1开启软件平台MATLAB，开启MATLAB编辑窗口； 2根据各种数值解法步骤编写M文件3保存文件并运行；4观察运行结果(数值或图形)；5根据观察到的结果写出实验报告，并浅谈学习心得体会。四、实验要求与任务根据实验内容和步骤，完成以下具体实验，要求写出实验报告（实验目的问题数学模型算法与编程计算结果分析、检验和结论心得体会）1数据插值：（1）轮船的甲板成近似半椭圆面形为了得到甲板的面积，首先测量得到横向最大相间8.534米；然后等间距地测得纵向高度，自左向右分别为：0.914, 5.060, 7.772, 8.717, 9.083, 9.144, 9.083

23、, 8.992, 8.687, 7.376, 2.073，计算甲板的面积。(2) 山区地貌图 在某山区（平面区域（0，2800）（0，2400）内，单位：米）测得一些地点的高程（单位：米）如下表所示，试作出该山区的地貌图和等高线图。2400200016001200800400 01430 1450 1470 1320 1280 1200 1080 9401450 1480 1500 1550 1510 1430 1300 12001460 1500 1550 1600 1550 1600 1600 16001370 1500 1200 1100 1550 1600 1550 13801270

24、1500 1200 1100 1350 1450 1200 11501230 1390 1500 1500 1400 900 1100 10601180 1320 1450 1420 1400 1300 700 900Y/X 0 400 800 1200 1600 2000 2400 28002数据拟合Malthus人口指数增长模型中参数从17901980年间美国每隔10年的人口记录如下表：年 份1790180018101820183018401850人口(106)3.95.37.29.612.917.123.2年 份1860187018801890190019101920人口(106)31.438.650.262.976.092.0106.5年 份193019401950196019701980 人口(106)123.2131.7150.7179.3204.0226.5 用以上数据检验马尔萨斯(Malthus)人口指数增长模型，根据检验结果进一步讨论马尔萨斯人口模型的改进。提示：Malthus 模型的基本假设是：人口的增长率为常数，记为 r。记时刻t的人口为x(t），（即x(t）为模型的状态变量）且初始时刻的人口为x0，于是得到如下微分方程：需要先求微分方程的解，再用数据拟合模型中的参数。