matlab的数据拟合与插值

matlab的数据拟合与插值
Matlab 的数据的分析处理-拟合与插值
在数学建模过程中，常常需要确定⼀个变量依存于另⼀个或更多的变量的关系，即确定这些变量之间的函数关系。

但在实际中确定这些变量之间函数函数关系时往往没有先验的依据，只能在收集的实际数据的基础上对若⼲合乎理论的形式进⾏试验，从中选择⼀个最有可能反映实际的函数形式，这就是统计学中的拟合和回归⽅程问题。

本节我们主要介绍如何分析处理实际中得到的数据。

下⾯先看⼀个例⼦。

例1 “⼈⼝问题”是我国最⼤社会问题之⼀，估计⼈⼝数量和发展趋势是我们制定⼀系列相关政策的基础。

有⼈⼝统计年鉴，可查到我国从1949年⾄1994
⼀般地，我们采⽤下⾯的分析处理⽅法：⾸先，在直⾓坐标系上作出⼈⼝数与年份的散点图象。

观察随着年份的增加⼈⼝数与年份变化关系，初步估计出他们之间的关系可近似地可看做⼀条直线。

那么我们如何把这条直线⽅程确定出来呢？并⽤他来估计1999年我国的⼈⼝数。

⽅法⼀：先选择能反映直线变化的两个点，如（1949，541.67），（1984，1034.75）⼆点确定⼀条直线，⽅程为 N = 14.088 t – 26915.842 ，代⼊t =1999,得N ≈12.46亿
⽅法⼆：可以多取⼏组点对，确定⼏条直线⽅程，将t = 1999代⼊，分别求出⼈⼝数，在取其算数平值。

⽅法三：可采⽤“最⼩⼆乘法”求出直线⽅程。

这就是曲线拟合的问题。

⽅法⼀与⽅法⼆都具有⼀定的局限性，下⾯我们重点介绍数据的曲线拟合。

所谓曲线拟合是指给定平⾯上的n 个点(x i ,y i ),i=1,2,….,n,找出⼀条曲线使之与这些点相当吻合，这个过程称之为曲线拟合。

最常见的曲线拟合是使⽤多项式来作拟合曲线。

曲线拟合最常⽤的⽅法是最⼩⼆乘法。

其原理是求f(x),使
21
])([i n
i i y x f -=∑=δ达到最⼩。

matlab 提供了基本的多项式曲线拟合函数命令
polyfit
格式::polyfit(x,y,n)
说明:polyfit(x,y,n)是找n 次多项式p(x)的系数,这些系数满⾜在最⼩⼆乘法意义下p(x(i)) ~= y(i).
已知⼀组数据，⽤什么样的曲线拟合最好呢?可以根据散点图进⾏直观观察，在此基础上，选择⼏种曲线分别拟合，然后⽐较，观察那条曲线的最⼩⼆乘指标最⼩。

下⾯我们给出常⽤的曲线（下⾯的,x y 为变量，,a b 等为参数）直线：y ax b =+
多项式：121231....n n n n n y a x a x a x a x a ---=+++++ (⼀般情况下，n 不宜过⾼，n=2,3) 双曲线：y=a
y b x
=
+ 指数曲线：bx y ae = 幂函数：b y ax =
有些曲线的拟合，为了利⽤数学软件，在拟合前需作变量替换，化为对未知数的线性函数。

思考：如果根据经验，曲线是双曲线a
y b x
=+或指数曲线bx y ae =及幂函数b
y ax =等，如何利⽤matlab 的多项式拟合函数来作曲线拟合？
例2：在化学反应中，为研究某化合物的浓度随时间的变化规律。

测得⼀组数据
本题是⼀个可以⽤数据的曲线拟合来解决的问题。

下⾯是利⽤matlab 编的⼀段程序。

clear; ％录⼊数据 xy=[1 4 2 6.4 3 8.0 4 8.4 5 9.28 6 9.5 7 9.7 8 9.86 9 10
10 10.2 11 10.32 12 10.42 13 10.5 14 10.55 15 10.58 16 10.6]; x=xy(:,1); y=xy(:,2);
plot(x,y,'r*');％画出散点图，观察曲线⾛势
hold on;t=0:.3:10;pxdxs=polyfit(x,y,2);pxd=poly2str(pxdxs,'x')
pxdx=polyval(pxdxs,t);plot(t,pxdx,'-k') ⽅法2：解下述⽅程组：（这是超定⽅程组（⽅程个数⼤于未知数个数的⽅程），这个⽅程组没有普遍意义下的解，但可以在最⼩⼆乘法意义下求解）
++=++=++=++=++=++=++=++=++=++=++=++=++=++=++=++=2^16*16*6.102
^15*15*58.102^14*14*55.102^13*13*5.102^12*12*42.102^11*11*32.102^10*10*2.102^9*9*102^8*8*86.92^7*7*7.92^6*6*5.92^5*5*28.92^4*4*4.82^3*3*0.82^2*2*4.62^1*1*4c b a c b a c b a c b a c b a c b a c b a c b a c b a c b a c b a c b a c b a c b a c b a c b a
把它写成矩阵乘法的形式： y ＝a*[a,b,c]'
其中，a=[1,1,1;1 2 2^2;1 3 3^2;1 4 4^2;1 5 5^2;1 6 6^2;1 7 7^2; 1 8 8^2;1 9 9^2;1 10 10^2; …
1 11 11^2;1 1
2 12^2;1 1
3 13^2;1 1
4 14^2;1 1
5 15^2;1 1
6 16^2];
y=[4 6.4 8.0 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]';
于是，abc=a\y
例3: 给药⽅案⼀种新药⽤于临床之前，必须设计给药⽅案. 药物进⼊机体后⾎液输送到全⾝，在这个过程中不断地被吸收、分布、代谢，最终排出体外，药物在⾎液中的浓
度，即单位体积⾎液中的药物含量，称为⾎药浓度。

我们将整个机体看作⼀个房室，称为中⼼室，并且假设室内⾎药浓度是均匀的(这样建⽴的模型称为⼀室模型)。

快速静脉
注射后，中⼼室内浓度⽴即上升,然后迅速下降。

当浓度太低时，达不到预期的治疗效果,如果浓度太⾼，⼜可能导致药物中毒或副作⽤太强。

临床上，每种药物有⼀个最⼩有效浓度c 1和⼀个最⼤有效浓度c 2。

设计给药⽅案时，要使⾎液中药物浓度保持在c 1~c 2之间。

本问题设c 1=10，c 2=25(ug/ml). 要设计给药⽅案,必须知道给药后⾎药浓度随时
间变化的规律。

从实验和理论两⽅⾯着⼿.
现对某药物进⾏临床实验,对某⼈⽤快速静脉注射⽅式⼀次注⼊该药物300mg 后,
问题: 1. 在快速静脉注射的给药⽅式下，研究⾎药浓度（单位体积⾎液中的药物含量）的变化规律
2. 给定药物的最⼩有效浓度和最⼤治疗浓度，设计给药⽅案：每次注射剂量多⼤；间隔时间多长。

分析:通过数据分析可以看到,⾎药浓度与时间的变化关系，符合负指数变化规律
为了建⽴模型,
除了刚才我们所做的假设以外,我们再做如下假设: 药物排除速
率与⾎药浓度成正⽐，⽐例系数为(0)k k >, ⾎液容积为v 不发⽣改变, 当0t =时注射剂量d, ⾎药浓度⽴即达到/d v . 于是,我们可建⽴如下的⼀室模型:
(0)/dc
kt dt
c d v ?=-?
=? 从⽽()kt
d c t
e v
-=
本问题中,已知300d mg =,要寻求()c t 与t 的关系,需要求出常数k 和v ,这需要⽤到拟合.
下⾯⽤matlab 的多项式拟合函数来拟合出k 和v ,程序如下:(这⾥需要先对数据进⾏处理) 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)
v=d/exp(a(2))
计算结果:0.2347,15.02k v ==
设每次注射剂量d , 间隔时间τ,⾎药浓度c(t) 满⾜12()c c t c << 初次剂量0d 应加⼤.将
给药⽅案记为：{}0,,d d t ,那么给药⽅案应满⾜:
初始药物浓度不⼤于(最好是等于)2c ,下次给药前⾎液浓度不低于1c (取下界等于
).
于
是
02
2
,
()d v c d
v c c =
=-,12k c c e
τ-=,从⽽可解
得:0375.5,225.3, 3.9d d τ===
提出给药⽅案:初次注射375.5mg ，以后每四⼩时注射225.3mg 。

课后思考：
1，已知⼀组数据，如何判断⽤什么样的曲线拟合最好？可否直观化？
2.下表给出了15个各种年纪的⼈的⾝⾼和体重的数据，试建⽴模型给出⼈的⾝
⼆、多项式插值
在⼯程实践合科学试验中，常常需要从⼀组实验观测数据(x i ,y i ),i=1,2..,n 来揭⽰两个变量之间的关系，⼀般可以⽤⼀个近似的函数关系式y=f(x)来表⽰，这种关系式的产⽣有两种⽅法：⼀种是曲线拟合，⼀种是插值。

拟合的思想前⾯我们已经说过，插值则要求求得的函数在节点（即每⼀个观测点）处的函数值⼀定要满⾜y i =f(x i )，⽽拟合不⼀定有此要求，它考虑到观测数据受随机误差的影响，寻求整体误差最⼩。

matlab ⽤来作插值计算的函数有⼀维插值的interp1，⼆维插值的interp2，三维插值的interp3，……,下⾯我们主要介绍⼀维插值。

调⽤格式为：
yi=interp1(x,y,xi,’method ’)
其中x ，y 为插值点，yi 为在被插值点xi 处的插值结果，x 、y 为向量，method 表⽰插值⽅法，matlab 的插值⽅法有：’nearest ’ （最邻近插值）、’linear ’（线性插值）、
’spline ’（三次样条插值）、’cubic ’（⽴⽅插值），缺省时表⽰线条插值。

所有的插值⽅法都要求x 是单调的，并且xi 的取值不能超过x 的范围例3：在⼀天的24个⼩时内，从零点开始每间隔2⼩时测得的环境温度数据分别为（℃）：
12，9,9,1,0,18,27,25,20,18,15,13
推测下午1点的环境温度。

键⼊：
x=0:2:24;y=[12,9,9,10,18,24,28,27,25,20,18,15,13];x1=13;y1=interp1(x,y,x1,’spline ’)
输出的y1即为所求。

如要得到⼀天24⼩时的温度曲线，只需继续键⼊： xi=0:1/3600:24;yi=interp1(x,y,xi,’spline ’); plot(x,y,’o ’,xi,yi)
课后思考题：
1,⽤电压V =14伏的电池给电容器充电，电容器上t 时刻的电压满⾜：
)exp()()(0τt
V V V t v ---=，
其中0V 是电容器的初始电压，τ是充电常数。

试⽤下列数据确定0V 和τ。

你⽤的⽅法是，结果是0V = ，τ= 。

2,⽕车⾏驶的路程、速度数据如下所⽰，编程计算从静⽌开始20分钟内⾛过的
3，⼩型⽕箭初始质量为1200千克，其中包括900千克燃料。

⽕箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产⽣30000⽜顿的恒定推⼒。

当燃料⽤尽时引擎关闭。

设⽕箭上升的整个过程中，空⽓阻⼒与速度平⽅成正⽐，⽐例系数记作k 。

⽕箭升空过程的数学模型为
kx T mg t t x x =-+-≤≤== 其中)(t x 为⽕箭在时刻t 的⾼度，120015m t =-为⽕箭在时刻t 的质量，T
（=30000⽜顿）为推⼒，g (=9.8⽶/秒2)为重⼒加速度, t 1 (=900/15=60秒)为引擎关闭时刻。

今测得⼀组数据如下（t-时间（秒），x -⾼度（⽶），v -速度（⽶/秒））：
现有两种估计⽐例系数k 的⽅法：
1，⽤每⼀个数据(t,x,v )计算⼀个k 的估计值（共11个），再⽤它们来估计k 。

2，⽤这组数据拟合⼀个k 。

请你分别⽤这两种⽅法给出k的估计值，对⽅法进⾏评价，并且回答，能否认为空⽓阻⼒系数k=0.5（说明理由）。