数学建模与数学实验后勤工程学院数学教研室拟 合实验目的实验内容
2、掌握用数学软件求解拟合问题。
1、直观了解拟合基本内容。
1,拟合问题引例及基本理论。
4,实验作业。
2,用数学软件求解拟合问题。
3,应用实例拟 合
2.拟合的基本原理
1,拟合问题引例拟 合 问 题 引 例 1
温度 t(0C) 20.5 32.7 51.0 73.0 95.7
电阻 R(?) 765 826 873 942 1032
已知热敏电阻数据：
求 600C时的电阻 R。
20 40 60 80 100
700
800
900
1000
1100
设 R=at+b
a,b为待定系数拟 合 问 题 引 例 2
t (h) 0.25 0.5 1 1.5 2 3 4 6 8
c (?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)下的图形
0
0
()
,
ktc t c e
ck
为 待 定 系 数
0 2 4 6 8
10
0
10
1
10
2
MATLAB(aa1)
曲 线 拟 合 问 题 的 提 法已知一组（二维）数据，即平面上 n个点 （ xi,yi) i=1,…n,
寻求一个函数（曲线） y=f(x),使 f(x) 在某种准则下与所有数据点最为接近，即曲线拟合得最好。
+
+
+
+ +
+
+
+
+
x
y
y=f(x)
(xi,yi)
i
i 为点 （ xi,yi) 与 曲线 y=f(x) 的距离拟合与插值的关系说明：
函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上完全不同。
实例,下面数据是某次实验所得，希望得到 x和 f之间的关系？x 1 2 4 7 9 12 13 15 17
f 1,5 3,9 6,6 1 1,7 1 5,6 1 8,8 1 9,6 2 0,6 2 1,1
MATLAB(cn)
问题,给定一批数据点，需确定满足特定要求的曲线或曲面解决方案：
若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，就是 数据拟合,又称曲线拟合或曲面拟合。
若要求所求曲线（面）通过所给所有数据点，就是 插值问题 ；
最临近插值、线性插值、样条插值与曲线拟合结果：
0 2 4 6 8 10 12 14 16 18
0
5
10
15
20
25
òaêy?Yμ?
s p l i n e
èy′àê?2μ
0 2 4 6 8 10 12 14 16 18
0
5
10
15
20
25
òaêy?Yμ?
l i n e s t
èy′àê?2μ
024681012141618
0
5
10
15
20
25
òaêy?Yμ?
nearest èy′àê?2μ
曲线拟合问题最常用的解法 —— 线性最小二乘法的基本思路第一步,先选定一组函数 r1(x),r2(x),…r m(x),m<n,令
f(x)=a1r1(x)+a2r2(x)+ …+a mrm(x) （ 1）
其中 a1,a2,…a m 为待定系数。
第二步,确定 a1,a2,…a m 的准则（最小二乘准则）：
使 n个点 （ xi,yi) 与 曲线 y=f(x) 的距离?i 的平方和最小 。
记
)2(])([
])([),,(
2
1 1
2
1 1
2
21
iik
n
i
m
k
k
i
n
i
n
i
iim
yxra
yxfaaaJ
问题归结为，求 a1,a2,…a m 使 J(a1,a2,…a m) 最小。
线性最小二乘法的求解：预备知识超定方程组,方程个数大于未知量个数的方程组
)(
2211
11212111
mn
yararar
yararar
nmnmnn
mm
即 Ra=y
nmnmnn
m
y
y
y
a
a
a
rrr
rrr
R
11
21
11211
,,其中超定方程一般是不存在解的矛盾方程组。
如果有向量 a使得 达到最小，
则称 a为上述 超定方程的最小二乘解 。
2
1
2211 )( im
n
i
imii yararar
线性最小二乘法的求解定理,当 RTR可逆时，超定方程组（ 3）存在最小二乘解，
且即为方程组
RTRa=RTy ------正则（正规）方程组的解,a=(RTR)-1RTy
所以，曲线拟合的最小二乘法要解决的问题，实际上就是求以下超定方程组的最小二乘解的问题。
nmnmn
m
y
y
y
a
a
a
xrxr
xrxr
R
11
1
111
,,
)()(
)()(其中
Ra=y （ 3）
线性最小二乘拟合 f(x)=a1r1(x)+ …+a mrm(x)中函数 {r1(x),… rm(x)}的选取
1,通过机理分析建立数学模型来确定 f(x)；
+ +
+
+
+
+ +
+
+
+ + +
+ ++
+ +
+
+
+
+ +
+ ++
+
+
+ ++
f=a1+a2x f=a1+a2x+a3x2 f=a1+a2x+a3x2
f=a1+a2/x f=aebx f=ae-bx
2,将数据 (xi,yi) i=1,…n 作图，通过直观判断确定 f(x)：
用 MATLAB解拟合问题
1、线性最小二乘拟合
2、非线性最小二乘拟合用 MATLAB作线性最小二乘拟合
1,作多项式 f(x)=a1xm+ …+a mx+am+1拟合,可利用已有程序,
a=polyfit(x,y,m)
2.对超定方程组
)(11 nmyaR nmmn
可得最小二乘意义下的解。
，用 yRa \?
3.多项式在 x处的值 y的计算命令,y=polyval（ a,x）
输出拟合多项式系数
a=[a1,…,am,am+1]’ （ 数组）
输入同长度数组 X,Y
拟合多项式次数即要求 出二次多项式,
3221)( axaxaxf
中 的
1 2 3(,,)A a a a
使得,11
2
1
[ ( ) ] ii
i
f x y
最 小例 对下面一组数据作二次多项式拟合
x i 0,1 0,2 0,4 0,5 0,6 0,7 0,8 0,9 1
y i 1,9 7 8 3,2 8 6,1 6 7,3 4 7,6 6 9,5 8 9,4 8 9,3 0 1 1,2
2
11 1
12
222
2
32
1111 11
,,1
,,1
,,1
xx y
a
yxx
a
a
yxx
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=R\y'
1
1
11
2
11
1
2
1
xx
xx
R此时
MATLAB(zxec1)
解法 1,解超定方程的方法
2）计算结果,Ａ = [-9.8108,20.1293,-0.0317]
0317.01293.208108.9)( 2 xxxf
RA y
2）计算结果,Ａ = [-9.8108,20.1293,-0.0317]
解法 2,用多项式拟合的命令
MATLAB(zxec2)
0 0,2 0,4 0,6 0,8 1
-2
0
2
4
6
8
10
12
0317.01293.208108.9)( 2 xxxf
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') %作出数据点和拟合曲线的图形
1,lsqcurvefit
已知 数据点,xdata=（ xdata1,xdata2,…,xdatan）
ydata=（ ydata1,ydata2,…,ydatan）
用 MATLAB作非线性最小二乘拟合两个求非线性最小二乘拟合的函数：
lsqcurvefit,lsqnonlin。
相同点和不同点：两 个命令都要先建立 M-文件 fun.m，定义函数 f(x)，但定义 f(x)的方式不同，请参考例题 。
2
1
1 ( (,) )
2
n
ii
i
F x x d a t a y d a t a
最 小
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,lb,ub);
(3) x =lsqcurvefit (‘fun’,x0,xdata,ydata,lb,ub,options);
(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 2 212( ) ( ) ( ) ( ) ( )T nf x f x L f x f x f x L
2,lsqnonlin
已知数据点,xdata=（ xdata1,xdata2,…,xdatan）
ydata=（ ydata1,ydata2,…,ydatan）
输入格式：
1） x=lsqnonlin（‘ fun’,x0）；
2） x= lsqnonlin （‘ fun’,x0,lb,ub）；
3） x= lsqnonlin （‘ fun’,x0，,lb,ub,options）；
4） [x,options]= lsqnonlin （‘ fun’,x0,… ）；
5） [x,options,funval]= lsqnonlin （‘ fun’,x0,… ）；
说明,x= lsqnonlin（‘ fun’,x0,options）；
fun是一个事先建立的定义函数 f(x)的 M-文件，
自变量为 x 迭代初值选项见无约束优化
1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 6 0 0 7 0 0 8 0 0 9 0 0 1 0 0 0
4,5 4 4,9 9 5,3 5 5,6 5 5,9 0 6,1 0 6,2 6 6,3 9 6,5 0 6,5 9
jt
310jc
10
0,0 2 2
1
11m i n (,,) [ ]
22
jkt
j
j
F a b k a b e c?
例 2 用下面一组数据拟合中的参数 a,b,k
ktbeatc 2.0.0)(
该问题即解最优化问题：
MATLAB(fzxec1)
1）编写 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）输入命令
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];
x0=[0.2,0.05,0.05];
x=lsqcurvefit ('curvefun1',x0,tdata,cdata)
f= curvefun1(x,tdata)
F(x,tdata)=,x=(a,b,k)
Tktkt beabea ),,( 101 02.002.0
解法 1,用命令 lsqcurvefit
3） 运算结果,
f =0.0043 0.0051 0.0056 0.0059 0.0061
0.0062 0.0062 0.0063 0.0063 0.0063
x =0.0063 -0.0034 0.2542
4）结论,a=0.0063,b=-0.0034,k=0.2542
0,0 2 0,2 5 4 2
0,0 0 5 0 8 4
0,0 0 6 3 0,0 0 3 4
0,0 0 6 3 0,0 0 3 4
t
t
c t e
e
MATLAB(fzxec2)
1） 编写 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)- cdata
2）输入命令,
x0=[0.2,0.05,0.05];
x=lsqnonlin('curvefun2',x0)
f= curvefun2(x)
函数 curvefun2的自变量是 x，
cdata和 tdata是已知参数，故应将 cdata tdata的值写在
curvefun2.m中解法 2 用命令 lsqnonlin
101 0.020.02
1 1 0
,,
(,,)ktkt T
f x F x td a ta c ta d a
a b e c a b e c
x=（ a,b,k）
3）运算结果为
f =1.0e-003 *（ 0.2322 -0.1243 -0.2495 -0.2413
-0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792）
x =0.0063 -0.0034 0.2542
可以看出，两个命令的计算结果是相同的 。
4）结论,即拟合得 a=0.0063 b=-0.0034 k=0.2542
MATLAB解应用问题实例
1、电阻问题
2、给药方案问题
*3、水塔流量估计问题
MATLAB(dianzu1)
电阻问题温度 t(0C) 20.5 32.7 51.0 73.0 95.7
电阻 R(?) 765 826 873 942 1032
例,由数据 拟合 R=a1t+a2
方法 1.用命令 polyfit(x,y,m)
得到 a1=3.3940,a2=702.4918
方法 2.直接用
yRa \?
结果相同。
MATLAB(dianzu2)
一室模型,将整个机体看作一个房室，称 中心室,
室内血药浓度是均匀的。快速静脉注射后，浓度立即上升；然后迅速下降。当浓度太低时，达不到预期的治疗效果；当浓度太高，又可能导致药物中毒或副作用太强。临床上，每种药物有一个最小有效浓度 c1和一个最大有效浓度 c2。设计给药方案时，要使血药浓度 保持在 c1~c2之间。本题设 c1=10，
c2=25(ug/ml).
拟 合 问 题 实 例 2给药方案 ——
一种新药用于临床之前，必须设计给药方案,
药物进入机体后血液输送到全身，在这个过程中不断地被吸收、分布、代谢，最终排出体外，药物在血液中的浓度，即单位体积血液中的药物含量，
称为 血药浓度。
在实验方面,对某人用快速静脉注射方式一次注入该药物 300mg后,在一定时刻 t(小时 )采集血药,测得血药浓度 c(ug/ml)如下表,
t (h) 0.25 0.5 1 1.5 2 3 4 6 8
c (?g/ml) 19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01
要设计给药方案,必须知道给药后血药浓度随时间变化的规律。从实验和理论两方面着手：
给药方案
1,在快速静脉注射的给药方式下，研究血药浓度（单位体积血液中的药物含量）的变化规律。
t
c2
c
c1
0?
问题
2,给定药物的最小有效浓度和最大治疗浓度，设计给药方案：每次注射剂量多大；间隔时间多长。
分析
理论：用一室模型研究血药浓度变化规律
实验：对血药浓度数据作拟合，符合负指数变化规律
3.血液容积 v,t=0注射剂量 d,血药浓度立即为 d/v.
2.药物排除速率与血药浓度成正比，比例系数 k(>0)
模型假设
1,机体看作一个房室，室内血药浓度均匀 —— 一室模型模型建立
d/c ( 0 ) 3?得：由假设
- k c
dt
dc 2?得：由假设
kte
v
dtc
)(
在此,d=300mg,t及 c（ t）在某些点处的值见前表，
需经拟合求出参数 k,v
用线性最小二乘拟合 c(t)
kte
v
dtc)(
)/l n (,,ln 21 vdakacy
ktvdc )/l n (ln
2/,
1
21
aedvak
atay
MATLAB(lihe1)
计算结果：
)(02.15),/1(2 3 4 7.0 lvhk
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))
程序：
用非线性最小二乘拟合 c(t)
给药方案 设计 cc2
c1
0? t
设每次注射剂量 D,间隔时间
血药浓度 c(t) 应 c1? c(t)? c2
初次剂量 D0 应加大
,,0 DD给药方案记为：
kecc 21
1
2ln1
c
c
k
2、
)(,1220 ccDcD
1、
计算结果,9.3,3.225,5.375
0DD
)(4),(225),(3750 hmgDmgD
给药方案：
c1=10,c2=25
k=0.2347
v=15.02
故可制定给药方案：
)(4),(225),(3750 hmgDmgD
即,
首次注射 375mg，
其余每次注射 225mg，
注射的间隔时间为 4小时。
估计水塔的流量
2,解题思路
3,算法设计与编程
1,问题某居民区有一供居民用水的园柱形水塔,一般可以通过测量其水位来估计水的流量,但面临的困难是,当水塔水位下降到设定的最低水位时,
水泵自动启动向水塔供水,到设定的最高水位时停止供水,这段时间无法测量水塔的水位和水泵的供水量,通常水泵每天供水一两次,每次约两小时,
水塔是一个高 12.2米,直径 17.4米的正园柱,按照设计,水塔水位降至约 8.2米时,水泵自动启动,水位升到约 10.8米时水泵停止工作,
表 1 是某一天的水位测量记录,试估计任何时刻
（ 包括水泵正供水时 ） 从水塔流出的水流量,及一天的总用水量,
表 1 水位测量记录
（符号 // 表示水泵启动）
时刻 ( h)
水位 ( c m )
0 0,92 1,84 2,95 3,87 4,98 5,90 7,01 7,93 8,97
968 948 931 913 898 881 869 852 839 822
时刻 ( h)
水位 ( c m )
9,98 10,92 10,95 12,03 12,95 13,88 14,98 15,90 16,83 17,93
/ / / / 1082 1050 1021 994 965 941 918 892
时刻 ( h)
水位 ( c m )
19,04 19,96 20,84 22,01 22,96 23,88 24,99 25,91
866 843 822 / / / / 1059 1035 1018
流量估计的解题思路拟合水位 ~时间函数确定流量 ~时间函数 估计一天总用水量拟合水位 ~时间函数测量记录看,一天有两个供水时段 （ 以下称第 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时刻的流量
MATLAB(llgj1)
4） 流量函数为,1 0 7 9.227 1 7 3.22 3 5 6.0)( 2 tttf
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时刻的流量
MATLAB(llgj2)
4） 流量函数为,8 3 1 3.17 5 1 2.87 5 2 9.00 1 8 6.0)( 23 ttttf
3,拟合供水时段的流量在第 1供水时段 （ t=9~11） 之前 （ 即第 1时段 ） 和之后 （ 即第 2时段 ） 各取几点,其流量已经得到,用它们拟合第 1供水时段的流量,为使流量函数在 t=9和 t=11连续,我们简单地只取 4个点,拟合
3次多项式 （ 即曲线必过这 4个点 ）,实现如下：
xx1=-polyval（ a1,[8 9]） ； %取第 1时段在 t=8,9的流量
xx2=-polyval（ a2,[11 12]） ； %取第 2时段在 t=11,12的流量
xx12=[xx1 xx2]；
c12=polyfit（ [8 9 11 12],xx12,3） ； %拟合 3次多项式
tp12=9,0.1,11；
x12=polyval（ c12,tp12） ； % x12输出第 1供水时段各时刻的流量
MATLAB(llgj3)
拟合的流量函数为：
078.3555879.737207.3)( 2 tttf
在第 2供水时段之前取 t=20,20.8两点的流水量，在该时刻之后（第 3时段）仅有 3个水位记录，我们用差分得到流量，然后用这 4个数值拟合第 2供水时段的流量如下：
dt3=diff（ t(22,24)）； %最后 3个时刻的两两之差
dh3=diff（ h(22,24)）； %最后 3个水位的两两之差
dht3=-dh3./dt3； %t(22)和 t(23)的流量
t3=[20 20.8 t(22) t(23)]；
xx3=[-polyval(a2,t3(1,2),dht3)]； %取 t3各时刻的流量
c3=polyfit（ t3,xx3,3）； %拟合 3次多项式
t3=20.8,0.1,24；
x3=polyval（ c3,tp3）； % x3输出第 2供水时段
（外推至 t=24）各时刻的流量
MATLAB(llgj4)
拟合的流量函数为,8 2 8 3.913 0 7 7.71 4 0 5.0)( 2 tttf
3,一天总用水量的估计第 1,2时段和第 1,2供水时段流量的积分之和,就是一天总用水量,虽然诸时段的流量已表为多项式函数,积分可以解析地算出,
这里仍用数值积分计算如下：
y1=0.1*trapz(x1)； %第 1时段用水量 （ 仍按高度计 ）,0.1为积分步长
y2=0.1*trapz(x2)； %第 2时段用水量
y12=0.1*trapz(x12)； %第 1供水时段用水量
y3=0.1*trapz(x3)； %第 2供水时段用水量
y=（ y1+y2+y12+y3） *237.8*0.01； %一天总用水量 （ ）
计算结果,y1=146.2,y2=266.8,y12=47.4,y3=77.3,y=1250.4
Lm 33 10?
MATLAB(llgjz)
4,流量及总用水量的检验计算出的 各时刻的流量 可用水位记录的数值微分来检验,用水量 y1可用第 1时段水位测量记录中下降高度 968-822=146来检验,类似地,y2用 1082-822=260检验,
供水时段流量 的一种 检验方法 如下：供水时段的用水量加上水位上升值 260是该时段泵入的水量,除以时段长度得到水泵的功率 （ 单位时间泵入的水量 ）,而两个供水时段水泵的功率应大致相等,第 1,2时段水泵的功率可计算如下：
p1=(y12+260)/2； %第 1供水时段水泵的功率
（ 水量仍以高度计 ）
tp4=20.8,0.1,23；
xp2=polyval（ c3,tp4） ； % xp2输出第 2供水时段各时刻的流量
p2=(0.1*trapz(xp2)+260)/2.2； %第 2供水时段水泵的功率
（ 水量仍以高度计 ）
计算结果,p1=154.5,p2=140.1 MATLAB (ll)
计算结果
( n 1,n 2 ) y 1,y 2,y 1 2,y 3 y p 1 p 2
( 3,4 ) 1 4 6,2,2 6 6,8,4 7,4,7 7,3 1 2 5 0,4 1 5 4,5 1 4 0,1
( 5,6 ) 1 4 6,5,2 5 7,8,4 6,1,7 6,3 1 2 8 2,4 1 5 3,7 1 4 0,1
2421 8 2 8 3.91 3 0 7 7.7 1 4 0 5.0
2111 8 3 1 3.1 7 5 1 2.8 7 5 2 9.0 0 1 8 6.0
119 0 7 8.3 5 5 5 8 7 9.73 7 2 0 7.3
90 1 0 7 9.22 7 1 7 3.2 2 3 5 6.0
)(
2
23
2
2
ttt
tttt
ttt
ttt
tf
流量函数为：
流量曲线见图
0 5 10 15 20 25
12
14
16
18
20
22
24
26
28
30
32
34
hour
cm
/h
o
u
o
r
0 5 10 15 20 25
12
14
16
18
20
22
24
26
28
30
32
34
hour
cm
/h
o
u
o
r
n=(3,4) n=(5,6)
练习 1 用给定的多项式，如 y=x3-6x2+5x-3，
产生一组数据 (xi,yi,i=1,2,…,n),再在 yi上添加随机干扰 (可用 rand产生 (0,1)均匀分布随机数,或用 rands产生 N(0,1)分布随机数 )，然后用 xi和添加了随机干扰的 yi作的 3次多项式拟合，与原系数比较。
如果作 2或 4次多项式拟合,结果如何？
练习 2,用电压 V=10伏的电池给电容器充电，
电容器上 t时刻的电压为,其中
V0是电容器的初始电压,是充电常数。试由下面一组 t,V数据确定 V0,。
t
eVVVtv )()( 0
t ( 秒 ) 0,5 1 2 3 4 5 7 9
V ( 伏 ) 6,36 6,48 7,26 8,22 8,66 8,99 9,43 9,63
用非线性最小二乘拟合 c(t)-用 lsqcurvefit
2、主程序 lihe2.m如下
clear
tdata=[0.25 0.5 1 1.5 2 3 4 6 8];
cdata=[19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01];
x0=[10,0.5];
x=lsqcurvefit('curvefun3',x0,tdata,cdata);
f=curvefun3(x,tdata)
x
MATLAB(lihe2)
1、用 M-文件 curvefun3.m定义函数
function f=curvefun3(x,tdata)
d=300
f=(x(1)\d)*exp(-x(2)*tdata)
% x(1)=v; x(2)=k
kte
v
dtc)(
2、掌握用数学软件求解拟合问题。
1、直观了解拟合基本内容。
1,拟合问题引例及基本理论。
4,实验作业。
2,用数学软件求解拟合问题。
3,应用实例拟 合
2.拟合的基本原理
1,拟合问题引例拟 合 问 题 引 例 1
温度 t(0C) 20.5 32.7 51.0 73.0 95.7
电阻 R(?) 765 826 873 942 1032
已知热敏电阻数据：
求 600C时的电阻 R。
20 40 60 80 100
700
800
900
1000
1100
设 R=at+b
a,b为待定系数拟 合 问 题 引 例 2
t (h) 0.25 0.5 1 1.5 2 3 4 6 8
c (?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)下的图形
0
0
()
,
ktc t c e
ck
为 待 定 系 数
0 2 4 6 8
10
0
10
1
10
2
MATLAB(aa1)
曲 线 拟 合 问 题 的 提 法已知一组（二维）数据，即平面上 n个点 （ xi,yi) i=1,…n,
寻求一个函数（曲线） y=f(x),使 f(x) 在某种准则下与所有数据点最为接近，即曲线拟合得最好。
+
+
+
+ +
+
+
+
+
x
y
y=f(x)
(xi,yi)
i
i 为点 （ xi,yi) 与 曲线 y=f(x) 的距离拟合与插值的关系说明：
函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上完全不同。
实例,下面数据是某次实验所得，希望得到 x和 f之间的关系？x 1 2 4 7 9 12 13 15 17
f 1,5 3,9 6,6 1 1,7 1 5,6 1 8,8 1 9,6 2 0,6 2 1,1
MATLAB(cn)
问题,给定一批数据点，需确定满足特定要求的曲线或曲面解决方案：
若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，就是 数据拟合,又称曲线拟合或曲面拟合。
若要求所求曲线（面）通过所给所有数据点，就是 插值问题 ；
最临近插值、线性插值、样条插值与曲线拟合结果：
0 2 4 6 8 10 12 14 16 18
0
5
10
15
20
25
òaêy?Yμ?
s p l i n e
èy′àê?2μ
0 2 4 6 8 10 12 14 16 18
0
5
10
15
20
25
òaêy?Yμ?
l i n e s t
èy′àê?2μ
024681012141618
0
5
10
15
20
25
òaêy?Yμ?
nearest èy′àê?2μ
曲线拟合问题最常用的解法 —— 线性最小二乘法的基本思路第一步,先选定一组函数 r1(x),r2(x),…r m(x),m<n,令
f(x)=a1r1(x)+a2r2(x)+ …+a mrm(x) （ 1）
其中 a1,a2,…a m 为待定系数。
第二步,确定 a1,a2,…a m 的准则（最小二乘准则）：
使 n个点 （ xi,yi) 与 曲线 y=f(x) 的距离?i 的平方和最小 。
记
)2(])([
])([),,(
2
1 1
2
1 1
2
21
iik
n
i
m
k
k
i
n
i
n
i
iim
yxra
yxfaaaJ
问题归结为，求 a1,a2,…a m 使 J(a1,a2,…a m) 最小。
线性最小二乘法的求解：预备知识超定方程组,方程个数大于未知量个数的方程组
)(
2211
11212111
mn
yararar
yararar
nmnmnn
mm
即 Ra=y
nmnmnn
m
y
y
y
a
a
a
rrr
rrr
R
11
21
11211
,,其中超定方程一般是不存在解的矛盾方程组。
如果有向量 a使得 达到最小，
则称 a为上述 超定方程的最小二乘解 。
2
1
2211 )( im
n
i
imii yararar
线性最小二乘法的求解定理,当 RTR可逆时，超定方程组（ 3）存在最小二乘解，
且即为方程组
RTRa=RTy ------正则（正规）方程组的解,a=(RTR)-1RTy
所以，曲线拟合的最小二乘法要解决的问题，实际上就是求以下超定方程组的最小二乘解的问题。
nmnmn
m
y
y
y
a
a
a
xrxr
xrxr
R
11
1
111
,,
)()(
)()(其中
Ra=y （ 3）
线性最小二乘拟合 f(x)=a1r1(x)+ …+a mrm(x)中函数 {r1(x),… rm(x)}的选取
1,通过机理分析建立数学模型来确定 f(x)；
+ +
+
+
+
+ +
+
+
+ + +
+ ++
+ +
+
+
+
+ +
+ ++
+
+
+ ++
f=a1+a2x f=a1+a2x+a3x2 f=a1+a2x+a3x2
f=a1+a2/x f=aebx f=ae-bx
2,将数据 (xi,yi) i=1,…n 作图，通过直观判断确定 f(x)：
用 MATLAB解拟合问题
1、线性最小二乘拟合
2、非线性最小二乘拟合用 MATLAB作线性最小二乘拟合
1,作多项式 f(x)=a1xm+ …+a mx+am+1拟合,可利用已有程序,
a=polyfit(x,y,m)
2.对超定方程组
)(11 nmyaR nmmn
可得最小二乘意义下的解。
，用 yRa \?
3.多项式在 x处的值 y的计算命令,y=polyval（ a,x）
输出拟合多项式系数
a=[a1,…,am,am+1]’ （ 数组）
输入同长度数组 X,Y
拟合多项式次数即要求 出二次多项式,
3221)( axaxaxf
中 的
1 2 3(,,)A a a a
使得,11
2
1
[ ( ) ] ii
i
f x y
最 小例 对下面一组数据作二次多项式拟合
x i 0,1 0,2 0,4 0,5 0,6 0,7 0,8 0,9 1
y i 1,9 7 8 3,2 8 6,1 6 7,3 4 7,6 6 9,5 8 9,4 8 9,3 0 1 1,2
2
11 1
12
222
2
32
1111 11
,,1
,,1
,,1
xx y
a
yxx
a
a
yxx
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=R\y'
1
1
11
2
11
1
2
1
xx
xx
R此时
MATLAB(zxec1)
解法 1,解超定方程的方法
2）计算结果,Ａ = [-9.8108,20.1293,-0.0317]
0317.01293.208108.9)( 2 xxxf
RA y
2）计算结果,Ａ = [-9.8108,20.1293,-0.0317]
解法 2,用多项式拟合的命令
MATLAB(zxec2)
0 0,2 0,4 0,6 0,8 1
-2
0
2
4
6
8
10
12
0317.01293.208108.9)( 2 xxxf
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') %作出数据点和拟合曲线的图形
1,lsqcurvefit
已知 数据点,xdata=（ xdata1,xdata2,…,xdatan）
ydata=（ ydata1,ydata2,…,ydatan）
用 MATLAB作非线性最小二乘拟合两个求非线性最小二乘拟合的函数：
lsqcurvefit,lsqnonlin。
相同点和不同点：两 个命令都要先建立 M-文件 fun.m，定义函数 f(x)，但定义 f(x)的方式不同，请参考例题 。
2
1
1 ( (,) )
2
n
ii
i
F x x d a t a y d a t a
最 小
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,lb,ub);
(3) x =lsqcurvefit (‘fun’,x0,xdata,ydata,lb,ub,options);
(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 2 212( ) ( ) ( ) ( ) ( )T nf x f x L f x f x f x L
2,lsqnonlin
已知数据点,xdata=（ xdata1,xdata2,…,xdatan）
ydata=（ ydata1,ydata2,…,ydatan）
输入格式：
1） x=lsqnonlin（‘ fun’,x0）；
2） x= lsqnonlin （‘ fun’,x0,lb,ub）；
3） x= lsqnonlin （‘ fun’,x0，,lb,ub,options）；
4） [x,options]= lsqnonlin （‘ fun’,x0,… ）；
5） [x,options,funval]= lsqnonlin （‘ fun’,x0,… ）；
说明,x= lsqnonlin（‘ fun’,x0,options）；
fun是一个事先建立的定义函数 f(x)的 M-文件，
自变量为 x 迭代初值选项见无约束优化
1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 6 0 0 7 0 0 8 0 0 9 0 0 1 0 0 0
4,5 4 4,9 9 5,3 5 5,6 5 5,9 0 6,1 0 6,2 6 6,3 9 6,5 0 6,5 9
jt
310jc
10
0,0 2 2
1
11m i n (,,) [ ]
22
jkt
j
j
F a b k a b e c?
例 2 用下面一组数据拟合中的参数 a,b,k
ktbeatc 2.0.0)(
该问题即解最优化问题：
MATLAB(fzxec1)
1）编写 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）输入命令
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];
x0=[0.2,0.05,0.05];
x=lsqcurvefit ('curvefun1',x0,tdata,cdata)
f= curvefun1(x,tdata)
F(x,tdata)=,x=(a,b,k)
Tktkt beabea ),,( 101 02.002.0
解法 1,用命令 lsqcurvefit
3） 运算结果,
f =0.0043 0.0051 0.0056 0.0059 0.0061
0.0062 0.0062 0.0063 0.0063 0.0063
x =0.0063 -0.0034 0.2542
4）结论,a=0.0063,b=-0.0034,k=0.2542
0,0 2 0,2 5 4 2
0,0 0 5 0 8 4
0,0 0 6 3 0,0 0 3 4
0,0 0 6 3 0,0 0 3 4
t
t
c t e
e
MATLAB(fzxec2)
1） 编写 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)- cdata
2）输入命令,
x0=[0.2,0.05,0.05];
x=lsqnonlin('curvefun2',x0)
f= curvefun2(x)
函数 curvefun2的自变量是 x，
cdata和 tdata是已知参数，故应将 cdata tdata的值写在
curvefun2.m中解法 2 用命令 lsqnonlin
101 0.020.02
1 1 0
,,
(,,)ktkt T
f x F x td a ta c ta d a
a b e c a b e c
x=（ a,b,k）
3）运算结果为
f =1.0e-003 *（ 0.2322 -0.1243 -0.2495 -0.2413
-0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792）
x =0.0063 -0.0034 0.2542
可以看出，两个命令的计算结果是相同的 。
4）结论,即拟合得 a=0.0063 b=-0.0034 k=0.2542
MATLAB解应用问题实例
1、电阻问题
2、给药方案问题
*3、水塔流量估计问题
MATLAB(dianzu1)
电阻问题温度 t(0C) 20.5 32.7 51.0 73.0 95.7
电阻 R(?) 765 826 873 942 1032
例,由数据 拟合 R=a1t+a2
方法 1.用命令 polyfit(x,y,m)
得到 a1=3.3940,a2=702.4918
方法 2.直接用
yRa \?
结果相同。
MATLAB(dianzu2)
一室模型,将整个机体看作一个房室，称 中心室,
室内血药浓度是均匀的。快速静脉注射后，浓度立即上升；然后迅速下降。当浓度太低时，达不到预期的治疗效果；当浓度太高，又可能导致药物中毒或副作用太强。临床上，每种药物有一个最小有效浓度 c1和一个最大有效浓度 c2。设计给药方案时，要使血药浓度 保持在 c1~c2之间。本题设 c1=10，
c2=25(ug/ml).
拟 合 问 题 实 例 2给药方案 ——
一种新药用于临床之前，必须设计给药方案,
药物进入机体后血液输送到全身，在这个过程中不断地被吸收、分布、代谢，最终排出体外，药物在血液中的浓度，即单位体积血液中的药物含量，
称为 血药浓度。
在实验方面,对某人用快速静脉注射方式一次注入该药物 300mg后,在一定时刻 t(小时 )采集血药,测得血药浓度 c(ug/ml)如下表,
t (h) 0.25 0.5 1 1.5 2 3 4 6 8
c (?g/ml) 19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01
要设计给药方案,必须知道给药后血药浓度随时间变化的规律。从实验和理论两方面着手：
给药方案
1,在快速静脉注射的给药方式下，研究血药浓度（单位体积血液中的药物含量）的变化规律。
t
c2
c
c1
0?
问题
2,给定药物的最小有效浓度和最大治疗浓度，设计给药方案：每次注射剂量多大；间隔时间多长。
分析
理论：用一室模型研究血药浓度变化规律
实验：对血药浓度数据作拟合，符合负指数变化规律
3.血液容积 v,t=0注射剂量 d,血药浓度立即为 d/v.
2.药物排除速率与血药浓度成正比，比例系数 k(>0)
模型假设
1,机体看作一个房室，室内血药浓度均匀 —— 一室模型模型建立
d/c ( 0 ) 3?得：由假设
- k c
dt
dc 2?得：由假设
kte
v
dtc
)(
在此,d=300mg,t及 c（ t）在某些点处的值见前表，
需经拟合求出参数 k,v
用线性最小二乘拟合 c(t)
kte
v
dtc)(
)/l n (,,ln 21 vdakacy
ktvdc )/l n (ln
2/,
1
21
aedvak
atay
MATLAB(lihe1)
计算结果：
)(02.15),/1(2 3 4 7.0 lvhk
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))
程序：
用非线性最小二乘拟合 c(t)
给药方案 设计 cc2
c1
0? t
设每次注射剂量 D,间隔时间
血药浓度 c(t) 应 c1? c(t)? c2
初次剂量 D0 应加大
,,0 DD给药方案记为：
kecc 21
1
2ln1
c
c
k
2、
)(,1220 ccDcD
1、
计算结果,9.3,3.225,5.375
0DD
)(4),(225),(3750 hmgDmgD
给药方案：
c1=10,c2=25
k=0.2347
v=15.02
故可制定给药方案：
)(4),(225),(3750 hmgDmgD
即,
首次注射 375mg，
其余每次注射 225mg，
注射的间隔时间为 4小时。
估计水塔的流量
2,解题思路
3,算法设计与编程
1,问题某居民区有一供居民用水的园柱形水塔,一般可以通过测量其水位来估计水的流量,但面临的困难是,当水塔水位下降到设定的最低水位时,
水泵自动启动向水塔供水,到设定的最高水位时停止供水,这段时间无法测量水塔的水位和水泵的供水量,通常水泵每天供水一两次,每次约两小时,
水塔是一个高 12.2米,直径 17.4米的正园柱,按照设计,水塔水位降至约 8.2米时,水泵自动启动,水位升到约 10.8米时水泵停止工作,
表 1 是某一天的水位测量记录,试估计任何时刻
（ 包括水泵正供水时 ） 从水塔流出的水流量,及一天的总用水量,
表 1 水位测量记录
（符号 // 表示水泵启动）
时刻 ( h)
水位 ( c m )
0 0,92 1,84 2,95 3,87 4,98 5,90 7,01 7,93 8,97
968 948 931 913 898 881 869 852 839 822
时刻 ( h)
水位 ( c m )
9,98 10,92 10,95 12,03 12,95 13,88 14,98 15,90 16,83 17,93
/ / / / 1082 1050 1021 994 965 941 918 892
时刻 ( h)
水位 ( c m )
19,04 19,96 20,84 22,01 22,96 23,88 24,99 25,91
866 843 822 / / / / 1059 1035 1018
流量估计的解题思路拟合水位 ~时间函数确定流量 ~时间函数 估计一天总用水量拟合水位 ~时间函数测量记录看,一天有两个供水时段 （ 以下称第 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时刻的流量
MATLAB(llgj1)
4） 流量函数为,1 0 7 9.227 1 7 3.22 3 5 6.0)( 2 tttf
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时刻的流量
MATLAB(llgj2)
4） 流量函数为,8 3 1 3.17 5 1 2.87 5 2 9.00 1 8 6.0)( 23 ttttf
3,拟合供水时段的流量在第 1供水时段 （ t=9~11） 之前 （ 即第 1时段 ） 和之后 （ 即第 2时段 ） 各取几点,其流量已经得到,用它们拟合第 1供水时段的流量,为使流量函数在 t=9和 t=11连续,我们简单地只取 4个点,拟合
3次多项式 （ 即曲线必过这 4个点 ）,实现如下：
xx1=-polyval（ a1,[8 9]） ； %取第 1时段在 t=8,9的流量
xx2=-polyval（ a2,[11 12]） ； %取第 2时段在 t=11,12的流量
xx12=[xx1 xx2]；
c12=polyfit（ [8 9 11 12],xx12,3） ； %拟合 3次多项式
tp12=9,0.1,11；
x12=polyval（ c12,tp12） ； % x12输出第 1供水时段各时刻的流量
MATLAB(llgj3)
拟合的流量函数为：
078.3555879.737207.3)( 2 tttf
在第 2供水时段之前取 t=20,20.8两点的流水量，在该时刻之后（第 3时段）仅有 3个水位记录，我们用差分得到流量，然后用这 4个数值拟合第 2供水时段的流量如下：
dt3=diff（ t(22,24)）； %最后 3个时刻的两两之差
dh3=diff（ h(22,24)）； %最后 3个水位的两两之差
dht3=-dh3./dt3； %t(22)和 t(23)的流量
t3=[20 20.8 t(22) t(23)]；
xx3=[-polyval(a2,t3(1,2),dht3)]； %取 t3各时刻的流量
c3=polyfit（ t3,xx3,3）； %拟合 3次多项式
t3=20.8,0.1,24；
x3=polyval（ c3,tp3）； % x3输出第 2供水时段
（外推至 t=24）各时刻的流量
MATLAB(llgj4)
拟合的流量函数为,8 2 8 3.913 0 7 7.71 4 0 5.0)( 2 tttf
3,一天总用水量的估计第 1,2时段和第 1,2供水时段流量的积分之和,就是一天总用水量,虽然诸时段的流量已表为多项式函数,积分可以解析地算出,
这里仍用数值积分计算如下：
y1=0.1*trapz(x1)； %第 1时段用水量 （ 仍按高度计 ）,0.1为积分步长
y2=0.1*trapz(x2)； %第 2时段用水量
y12=0.1*trapz(x12)； %第 1供水时段用水量
y3=0.1*trapz(x3)； %第 2供水时段用水量
y=（ y1+y2+y12+y3） *237.8*0.01； %一天总用水量 （ ）
计算结果,y1=146.2,y2=266.8,y12=47.4,y3=77.3,y=1250.4
Lm 33 10?
MATLAB(llgjz)
4,流量及总用水量的检验计算出的 各时刻的流量 可用水位记录的数值微分来检验,用水量 y1可用第 1时段水位测量记录中下降高度 968-822=146来检验,类似地,y2用 1082-822=260检验,
供水时段流量 的一种 检验方法 如下：供水时段的用水量加上水位上升值 260是该时段泵入的水量,除以时段长度得到水泵的功率 （ 单位时间泵入的水量 ）,而两个供水时段水泵的功率应大致相等,第 1,2时段水泵的功率可计算如下：
p1=(y12+260)/2； %第 1供水时段水泵的功率
（ 水量仍以高度计 ）
tp4=20.8,0.1,23；
xp2=polyval（ c3,tp4） ； % xp2输出第 2供水时段各时刻的流量
p2=(0.1*trapz(xp2)+260)/2.2； %第 2供水时段水泵的功率
（ 水量仍以高度计 ）
计算结果,p1=154.5,p2=140.1 MATLAB (ll)
计算结果
( n 1,n 2 ) y 1,y 2,y 1 2,y 3 y p 1 p 2
( 3,4 ) 1 4 6,2,2 6 6,8,4 7,4,7 7,3 1 2 5 0,4 1 5 4,5 1 4 0,1
( 5,6 ) 1 4 6,5,2 5 7,8,4 6,1,7 6,3 1 2 8 2,4 1 5 3,7 1 4 0,1
2421 8 2 8 3.91 3 0 7 7.7 1 4 0 5.0
2111 8 3 1 3.1 7 5 1 2.8 7 5 2 9.0 0 1 8 6.0
119 0 7 8.3 5 5 5 8 7 9.73 7 2 0 7.3
90 1 0 7 9.22 7 1 7 3.2 2 3 5 6.0
)(
2
23
2
2
ttt
tttt
ttt
ttt
tf
流量函数为：
流量曲线见图
0 5 10 15 20 25
12
14
16
18
20
22
24
26
28
30
32
34
hour
cm
/h
o
u
o
r
0 5 10 15 20 25
12
14
16
18
20
22
24
26
28
30
32
34
hour
cm
/h
o
u
o
r
n=(3,4) n=(5,6)
练习 1 用给定的多项式，如 y=x3-6x2+5x-3，
产生一组数据 (xi,yi,i=1,2,…,n),再在 yi上添加随机干扰 (可用 rand产生 (0,1)均匀分布随机数,或用 rands产生 N(0,1)分布随机数 )，然后用 xi和添加了随机干扰的 yi作的 3次多项式拟合，与原系数比较。
如果作 2或 4次多项式拟合,结果如何？
练习 2,用电压 V=10伏的电池给电容器充电，
电容器上 t时刻的电压为,其中
V0是电容器的初始电压,是充电常数。试由下面一组 t,V数据确定 V0,。
t
eVVVtv )()( 0
t ( 秒 ) 0,5 1 2 3 4 5 7 9
V ( 伏 ) 6,36 6,48 7,26 8,22 8,66 8,99 9,43 9,63
用非线性最小二乘拟合 c(t)-用 lsqcurvefit
2、主程序 lihe2.m如下
clear
tdata=[0.25 0.5 1 1.5 2 3 4 6 8];
cdata=[19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01];
x0=[10,0.5];
x=lsqcurvefit('curvefun3',x0,tdata,cdata);
f=curvefun3(x,tdata)
x
MATLAB(lihe2)
1、用 M-文件 curvefun3.m定义函数
function f=curvefun3(x,tdata)
d=300
f=(x(1)\d)*exp(-x(2)*tdata)
% x(1)=v; x(2)=k
kte
v
dtc)(
