温度分布的曲线拟合Word格式文档下载.docx
《温度分布的曲线拟合Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《温度分布的曲线拟合Word格式文档下载.docx（20页珍藏版）》请在冰点文库上搜索。
3.三次样条插值拟合
4.T7的三角多项式拟合
5.有4个控制点的贝塞尔曲线拟合
2.实验内容
、线性最小二乘拟合定理5.1（最小二乘拟合曲线）设｛（x^yk）｝：
#有N个点，其中横坐标｛兀｝仁是确定的
最小—乘拟合曲线
（1）
y=AxB
的系数是下列线性方程组的解，这些方程称为正规方程：
fN2）fN）N
lZXkIA+庄XkB=Ex<
yk*4丿丿7
NN
'
XkANB八yk2k/
核心代码为：
%求方程组am=b的根
m=a\b;
x1=1:
0.1:
24;
y1=m
（1）*x1+m
（2）;
%绘图，其中（x,y）为已知点，用红色的星号表示，y1为拟合曲线
plot（x,y,'
*r'
x1,y1）
gridon
legendf已知点'
'
最小二乘拟合'
）
主要算法为：
（1）.输入x,y；
NNNN
（2）.求正规方程的系数7x2?
Xk^yk?
Xkyk
k二k二k斗心
（3）.解正规方程组am=b（4）.绘制拟合曲线
图1线性的最小二乘拟合流程图
、曲线的最小二乘抛物线拟合
定理5.3（最小二乘抛物线拟合）设{（Xk,yQ}Nm有N个点，横坐标是确定的
二乘抛物线的系数表示为
y=f（x）二Ax2BxC
求解A,B和C的线性方程组为
『N\「N、（“）N
IZx：
IA+!
送X；
B+!
EX：
C=送ykx：
1心丿I心丿1心丿k¥
fN3）fN2）"
）N
.区xkIA+gxkB+：
瓦xkC=£
ykxk1心丿1心丿1心丿kT
NNN
|Hx：
A亠―xk|BNC八yk
kd?
kdkm
最小
⑶
⑷
根据式（4），核心代码为:
a（1,1）=sum（x.A4）;
a（2,3）=sum（x）;
b
（1）=（x.A2）*y'
;
b
（2）=x*y'
算法流程图为：
图2抛物线的最小二乘拟合流程图
三、三次样条插值拟合
定义5.1设｛（Xk,yj｝打有N-1个点，其中a=X。
：
为：
|1（:
:
Xn：
b。
如果存在N个三次多项式Sk（X），系数为Sk,o,Sk,i,Sk,2和Sk,3，满足如下性质：
IS（x）=Sk,0+Sk,i（x-xk）+Sk,2（X-Xk）+Sk,3（X-Xk）
II怒）二Yk
III.Sk（xk1^-Sk1（xk1）iv.Sk（xk.J=q1（xk1）
V.Sk'
（xk1）=Sk1（xk1）
x•[xk,xk1],k=0,11（,N-1k=O,1,|l|,N
（5）
k=O,1,|l|,N-2
k=0,1川，N—2
k=0,1,|l|,N-2
则称函数为三次样条函数。
令mk=s"
（Xk）,mk1=s"
（Xk1）,hk二x—-Xk和dk=山吐，可得包含口㈠皿和hk
mki的重要关系式：
hkjmkj2（hkj•hk）mkhEk1二山（6）
其中Uk=6（dk-dk」）,k=1,2」|（,N-1
方程组（6）中的未知数是要求的值｛叫｝，而且其他的项是可以通过数据点集｛（Xk.Yk）｝进行简单数学计算得到的常量。
因此方程组（6）是包含N1个未知数，具有N-1个线性方程组的不定方程组。
所以需要另外两个方程组才能求解。
可通过它们消去方程组（6）中的第一个方程的m0和第个方程的mN。
如果给定m。
，则可以计算出m°
ho，而且方程组⑹的第一个方程（当k=1时）为:
2（hoh）m1gm2二5-hom。
⑺
如果给定mN，则可以计算出hN^mN，而且方程组⑹的最后一个方程（当k=N-1时）为：
hN_2mN-2'
2（hNJ2hN4）mN4~UN4_hN-jmN（8）
考虑方程组⑹以及方程组⑺和方程组（8），其中k二2,3川I,N-2，可形成N-1阶线性方程组，包含系数m1,m2,,mN4。
重写方程组（6）中的方程1到方程N-1方程组HM=V，表示为：
b1
a1
当得到系数｛叫｝后，可以用如下公式计算Sk（x）的样条系数｛Sk,j｝
sk,o=yk，
Sk’^dk-h^g叫1）
sk,2-
mik
（10）
6hk
为了更有效地计算，每个三次多项式
Sk（x）可表示成嵌套形式：
（11）
Sk（x）=（（Sk,3WSk,2）wSk,i）wyk,其中w^x-
其中Sjx）在给区间Xk^x乞Xki内使用。
fork=2:
N-1
temp=A（k-1）/B（k-1）;
B（k）=B（k）-temp*C（k-1）;
U（k）=U（k）-temp*U（k-1）;
end
%求m（0）和m（N）
M
（1）=2*（D
（1）-dxO）/H
（1）-M
（2）/2;
M（N+1）=3*（dxn-D（N））/H（N）-M（N）/2;
%求样条系数s（k,j）
fork=0:
S（k+1,1）=（M（k+2）-M（k+1））/（6*H（k+1））;
S（k+1,2）=M（k+1）/2;
S（k+1,3）=D（k+1）-H（k+1）*（2*M（k+1）+M（k+2））/6;
S（k+1,4）=Y（k+1）;
图3三次样条拟合流程图
四、T7的三角多项式拟合
定义5.4具有如下形式的级数；
称为M阶的三角多项式。
定理5.8（离散傅里叶级数）设有N1个点｛（&
』"
｝：
£
，其中yj二f（x），而且横
坐标之间等距，即：
Xj=7M'
j,其中j=0,1」川，N（13）
Nn
如果f（x）的周期为2二，而且2M:
N，则存在式（12）所示的三角多项式Tm（x），
使得下式的值最小
N
、（f（Xk）-TM（x））
k4
多项式的系数
j和bj可通过如下公式计算：
2N
ajf（Xk）cos（jXk）,其中j=0,1,|）|,M
Nk4
bjf（Xk）sin（jXk）,其中j=0,1,|l（,M
（14）
（15）
（16）
图4三角多项式拟合流程图
核心代码为:
%计算A和B
forj=1:
M
A（j+1）=cos（j*X）*Y'
B（j+1）=sin（j*X）*Y'
%求三角多项式T
T=A
（1）;
T=T+A（j+1）*cos（j*x）+B（j+1）*sin（j*x）;
五、有4个控制点的贝塞尔曲线拟合
定义5.5N阶伯恩斯坦多项式定义为
N.B,N（t）二匚t，（1-t）亠（17）
廿+心N!
i=0,1,2,N,其中VUi!
（N-i）!
定义5.6给定一个控制点集，｛R｝寫，其中R，定义
P（t）八RBi,N（t）（18）
i=0
为N阶贝塞尔曲线，其中Bi,N（t）,i=0,1,,N，是N阶伯恩斯坦多项式，"
[0,1]。
公式（18）中的控制点是表示平面中x和y坐标的有序对。
可将控制点作为向量，儿对
应的伯恩斯坦多项式作为标量处理，
这样公式（18）可参数化表示为
P（t）=（x（t）,yt其中
X（t）二"
xiBi,N（t）
i=0
%B%3阶伯恩斯坦多项式系数矩阵
B=[（1-t）A3,3*t*（1-t）A2,3*tA2*（1-t）,tA3];
%xx为3阶贝塞尔曲线横坐标，yy为纵坐标xx=0;
yy=0;
fori=n+1:
n+4
xx=xx+x（i）*B（i-n）;
yy=yy+y（i）*B（i-n）;
图5贝塞尔曲线拟合流程图
3.实验结果及分析
从图6中看出由于温度变化快，数据比较分散，不宜用最小二乘直线进行拟合
图7最小二乘抛物线拟合曲线
比较图6和图7,显然图7中的拟合曲线更贴近于温度点坐标。
而图7能够直观的反映一天中的温度变化趋势。
但所拟合的曲线明显与数据点偏离较大。
从图8中可以很容易看出拟合曲线将经过样点，读出该地区一天的温度变化走势：
上午时段温度快速升高，到12点左右达到最高。
下午时段，温度慢慢降低，直到24时温度降至最低并保持一段时间。
从图9中可以很容易读出该地区一天的温度变化走势，与图8所反映的温度变化
规律大致相同。
与图三次样条曲线相比，三角多项式拟合出的曲线更光滑。
图10贝塞尔曲线拟合曲线
图8与图10所反映的温度变化规律大致相同，能够很容易地看出该地区一天温
度变化的走势。
图10采用4个控制点的贝塞尔拟合曲线在每段曲线交点出，其斜率变
化较大
4.结论
1.根据五种不同的拟合结果可知，由于线性最小二乘拟合和抛物线最小二乘拟合误差很大，三角多项式拟合次之，而三次样条拟合和贝塞尔曲线拟合最为精确，误差最小。
2.从本次实验可知，在选择最佳的拟合方式前可以先绘制出已知点的离散分布图像,根据其变化趋势合理选择最佳拟合方式。
3.在本实验中，三次样条拟合和贝塞尔曲线拟合最为精确，但两者相比，三次样条
拟合程序要复杂和繁琐些，而贝塞尔曲线拟合较为简单，故本实验中的最佳拟合方式为贝塞尔曲线拟合。
附件（代码）一、线性的最小二乘拟合
%时间x=1:
%温度y=[58,58,58,58,57,57,57,58,60,64,67,68,66,66,65,64,63,63,62,61,60,60,59,58];
%初始化正规方程的系数矩阵a,ba=zeros（2,2）;
b=zeros（2,1）;
a（1,1）=x*x'
a（1,2）=sum（x）;
a（2,1）=a（1,2）;
a（2,2）=24;
b
（1）=x*y'
b
（2）=sum（y）;
%求方程组am=b的根m=a\b;
二、曲线的最小二乘抛物线拟合；
clear
%时间
x=1:
%温度
y=[58,58,58,58,57,57,57,58,60,64,67,68,66,66,65,64,63,63,62,61,60,60,59,58];
%初始化正规方程的系数矩阵a,b
a=zeros（3,3）;
b=zeros（3,1）;
a（1,2）=sum（x.A3）;
a（1,3）=sum（x.A2）;
a（2,2）=a（1,3）;
a（2,3）=sum（x）;
a（3,1）=a（2,2）;
a（3,2）=a（2,3）;
a（3,3）=24;
b
（1）=（x.A2）*y'
b（3）=sum（y）;
y1=m
（1）*x142+m
（2）.*x1+m（3）;
%绘图，其中（x,y）为已知点，用红色的星号表示，y1最小二乘抛物线
最小二乘抛物线'
N=23;
X=1:
N+1;
Y=[58,58,58,58,57,57,57,58,60,64,67,68,66,66,65,64,63,63,62,61,60,60,59,58];
%dxO=S'
（xO）,dxn=S'
（xn）,为自由边界条件
dx0=0;
dxn=0;
H=diff（X）;
D=diff（Y）./H;
A=H（2:
N-1）;
B=2*（H（1:
N-1）+H（2:
N））;
C=H（2:
N）;
U=6*diff（D）;
%clampedsplineendpointconstraints
B
（1）=B
（1）-H
（1）/2;
U
（1）=U
（1）-3*（D
（1）-dx0）;
B（N-1）=B（N-1）-H（N）/2;
U（N-1）=U（N-1）-3*（dx0-D（N））;
M（N）=U（N-1）/B（N-1）;
fork=N-2:
-1:
M（k+1）=（U（k）-C（k）*M（k+2））/B（k）;
M⑴=2*（D⑴-dxO）/H
（1）-M
（2）/2;
%绘图，其中（X,Y）为已知点，用红色的星号表示，y为三次样条曲线
x=j:
0.01:
j+1;
y=polyval（S（j,:
）,x-X（j））;
plot（X,Y,'
x,y）,holdon
legend（'
已知点'
三次样条曲线'
symsx
M=7;
X=-pi:
2*pi/N:
pi;
%A为包含cos（jx）系数的矩阵
%B为包含sin（jx）系数的矩阵
A=zeros（1,M+1）;
B=zeros（1,M+1）;
%Yends为非连续点处的值
Yends=（丫⑴+Y（N））/2;
Y
（1）=Yends;
Y（N）=Yends;
A
（1）=sum（Y）;
A=2*A/N;
B=2*B/N;
A
（1）=A
（1）/2;
%绘图，其中（X,Y）为已知点，用红色的星号表示，y为
x1=-pi:
.01:
y1=subs（T,x,x1）;
x1,y1）
三角多项式曲线'
%时间，第25个值为构造值
24+1;
%温度，第25个值为构造值
y=[58,58,58,58,57,57,57,58,60,64,67,68,66,66,65,64,63,63,62,61,60,60,59,58,58];
symst
n=0;
fork=1:
8%xx为3阶贝塞尔曲线横坐标，yy为纵坐标
xx=0;
fori=n+1:
n=n+3;
tt=0:
1;
bx=subs（xx,t,tt）;
by=subs（yy,t,tt）;
%绘图，其中（x,y）为已知点，用红色的星号表示，by为贝塞尔曲线
bx,by）,holdon
贝塞尔曲线'
升级会员 