温度分布的曲线拟合

《数值计算方法》实验报告 1

温度分布的曲线拟合 学号:XX 姓名：XXX 1. 实验描述 美国洛杉矶郊区11月8日的温度（华氏温度）如表1所示。采用24小时制。 表1 温度数据 时间，p.m. 温度 时间，a.m. 温度 1 66 1 58 2 66 2 58 3 65 3 58 4 64 4 58 5 63 5 57 6 63 6 57 7 62 7 57 8 61 8 58 9 60 9 60 10 60 10 64 11 59 11 67 午夜 58 正午 68

要求：1.线性的最小二乘拟合 2.曲线的最小二乘抛物线拟合； 3.三次样条插值拟合 4.T7的三角多项式拟合 5.有4个控制点的贝塞尔曲线拟合 2. 实验内容

一、线性最小二乘拟合 定理5.1（最小二乘拟合曲线）设1{(,)}Nkkkxy有N个点，其中横坐标1{}Nkkx是确定的。《数值计算方法》实验报告 2

最小二乘拟合曲线 yAxB (1)

的系数是下列线性方程组的解，这些方程称为正规方程：

211111NNNkkkkkkkNNkkkkxAxBxyxANBy











 (2)

核心代码为： %求方程组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) grid on legend('已知点','最小二乘拟合')

主要算法为： (1).输入x,y； (2).求正规方程的系数21Nkkx,1Nkkx,1Nkky

,1Nkkkxy

(3).解正规方程组am=b (4).绘制拟合曲线 《数值计算方法》实验报告 3

二、曲线的最小二乘抛物线拟合 定理5.3（最小二乘抛物线拟合）设1{(,)}Nkkkxy有N个点，横坐标是确定的。最小二乘抛物线的系数表示为 2()yfxAxBxC

(3)

求解,AB和C的线性方程组为

432211113211112111NNNNkkkkkkkkkNNNNkkkkkkkkk

NNNkkkkkk

xAxBxCyxxAxBxCyxxAxBNCy (4)

根据式(4)，核心代码为： a(1,1)=sum(x.^4); a(2,3)=sum(x); b(1)=(x.^2)*y';

开始 输入x，y

21Nkkx,1Nkkx,1Nkky

,1Nkkkxy

解正规方程组am=b 绘图 结束 图1 线性的最小二乘拟合流程图 《数值计算方法》实验报告 4

b(2)=x*y'; %求方程组am=b的根 m=a\b; 算法流程图为：

三、三次样条插值拟合 定义5.1 设1{(,)}Nkkkxy有1N个点，其中01Naxxxb。如果存在N个三次多项式()kSx，系数为,0,1,2,,kkksss和,3ks，满足如下性质：

开始 输入x，y

21Nkkx,1Nkkx,1Nkky

,1Nkkkxy

解正规方程组am=b

绘图 结束

2()fxAxBxC

图2 抛物线的最小二乘拟合流程图 《数值计算方法》实验报告 5

23kk,0k,1kk,2kk,3k

1111''111''''111I.S(x)=s+s(x-x)+s(x-x)+s(x-x)[,],0,1,,1.()0,1,,.()()0,1,,2.()()0,1,,2.()()0,1,,2kkkkkkkkkkkkkkkk

xxxkNIISxykNIIISxSxkNIVSxSxkNVSxSxkN

(5)

则称函数为三次样条函数。 令''''11(),()kkkkmSxmSx，1kkkhxx和1kkkkyydh，可得包含1,kkmm和

1km的重要关系式：

11112()kkkkkkkkhmhhmhmu (6)

其中16(),1,2,,1kkkuddkN 方程组(6)中的未知数是要求的值{}km，而且其他的项是可以通过数据点集{(,)}kkxy进行简单数学计算得到的常量。因此方程组(6)是包含1N个未知数，具

有1N个线性方程组的不定方程组。所以需要另外两个方程组才能求解。可通过它们消去方程组(6)中的第一个方程的0m和第个方程的Nm。 如果给定0m，则可以计算出00mh，而且方程组(6)的第一个方程（当k=1时）为：

011121002()hhmhmuhm (7)

如果给定Nm，则可以计算出1NNhm，而且方程组(6)的最后一个方程（当k=N-1时）为：

22211112()NNNNNNNNhmhhmuhm (8)

考虑方程组(6)以及方程组(7)和方程组(8)，其中2,3,,2kN，可形成1N

阶线性方程组，包含系数121,,,Nmmm。 重写方程组(6)中的方程1到方程1N，得到一个包含121,,,Nmmm的三角线性方程组HMV，表示为：

1112222

32222211111NNNNNNNNN

bcmvabcmv

abcmvabmv

(9)

当得到系数{}km后，可以用如下公式计算()kSx的样条系数,{}kjS。 《数值计算方法》实验报告 6

1,0,11,2,3(2),6,26kkkkkk

kkkkkk

k

hkmmsysd

mmmssd

h

(10)

为了更有效地计算，每个三次多项式()kSx可表示成嵌套形式： ,3,2,1()(()),其中kkkkkkSxswswswywxx (11)

其中()kSx在给区间1kkxxx内使用。 核心代码为： for k=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)-dx0)/H(1)-M(2)/2; M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2; %求样条系数s(k,j) for k=0:N-1 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); end 算法流程图为： 《数值计算方法》实验报告 7

四、T7的三角多项式拟合 定义5.4 具有如下形式的级数； 0()(cos()sin())2MjjaTxajxbjx

(12)

称为M阶的三角多项式。 定理5.8 （离散傅里叶级数） 设有1N个点0{(,)}Nkkjxy，其中()jyfx，而且横坐标之间等距，即： 21,0,1,,其中jjxjN

Nn (13)

如果()fx的周期为2，而且2MN，则存在式(12)所示的三角多项式()MTx，使得下式的值最小。

开始 结束 绘图 1kkkkyydh 16()kkkudd



构造HM=V 样条系数,{}kjS

1kkkhxx 输入X，Y

图3 三次样条拟合流程图 《数值计算方法》实验报告 8

21(()())NkMkfxTx

(14)

多项式的系数ja和jb可通过如下公式计算：

12()cos(),0,1,,NjkkkafxjxjMN其中 (15)

12()sin(),0,1,,NjkkkbfxjxjMN其中 (16)

核心代码为： %计算A和B for j=1:M

开始 结束 输入X和Y

j>M A(j+1)<=cos(j*X)*Y' B(j+1)<=sin(j*X)*Y';

三角多项式T 绘图

j<=1 j<=j+1

图4 三角多项式拟合流程图 Y N