数值计算方法实验报告5—温度分布的曲线拟合

«数值计算方法»实验报告
标题：温度分布的曲线拟合
1.实验描述：
在科学技术工程和实验中，经常需要从大量的实验数据中寻找拟合曲线，最 简单的是一维情形（一元函数） ，此时数据的形式为 x 和 y 坐标的有序对，如： ( x1 , y1 ),..., ( xN , y N ), 这里的横坐标{ x}是明确的。 数值计算方法的目的之一是求解一个将自变量与因变量联系起来的拟合函 数。求解拟合函数的方法有多种，常见的方法有：线性最小二乘拟合、多项式拟 合（最小二乘抛物线拟合） 、样条插值拟合（三次样条拟合） 、三角多项式拟合、 贝塞尔曲线拟合这五种方法。 本次实验分别利用上述五种方法对一组温度数据进行拟合， 通过拟合的结果 比较这五种方法的优缺点（主要考虑误差） 。
aj =
2 2 f ( xk ) cos( jxk ), 其中; j = 0,1,..., M ; bj = ∑ f ( xk ) sin( jxk ), 其中:j = 1, 2,..., M ∑ N k =1 N k =1
a0 M + ∑ ( a j cos( jx ) + b j sin( jx )) 2 j =1
N
E ( A, B, C ) = ∑ ( Axk2 + Bx k + C − y k ) 2 , 分别对A, B, C求偏导数：
k =1
∂E ( A, B, C ) ∂E ( A, B, C ) ∂E ( A, B, C ) = 0， = 0, =0 A B C 最终得到上述正规方程。
③样条插值拟合（三次样条拟合）法：
N k =1 N k =1 N k =1 N k =1 N N
( ∑ xk2 ) A + ( ∑ xk ) B = ∑ xk yk
N k =1
( ∑ xk ) A + NB = ∑ yk
N N
解得： A =
∑ xk yk − N x y
k =1 N
∑(x
=
k =1 N
k
பைடு நூலகம்
− x)( yk − y) , B = y − Ax
关键代码 1（求解 A,B,C 及误差 E） ： D=[sum(X.^4) sum(X.^3) sum(X.^2);sum(X.^3) sum(X.^2) sum(X);sum(X.^2) sum(X) N]; F=[sum(T.*X.^2);sum(T.*X);sum(T)]; P=D\F A=P(1) B=P(2) C=P(3) E=sqrt(sum((A*X.^2+B*X+C-T).^2)/N) 关键代码 2（绘图） ： plot(X,A*X.^2+B*X+C) hold on plot(X,A*X.^2+B*X+C-T,'r') 拟合曲线及误差曲线如下： %画出拟合曲线 %在一幅图上画出拟合曲线及误差曲线 %画出误差曲线 %D 为正规方程的系数矩阵 %F 为正规方程的常数向量 %求解 A,B,C %A 的值 %B 的值 %C 的值 %计算均方根误差 E
设{( x , y )} 有N 个点，其中横坐标 {xk }k =1是确定的。最小二乘抛物线拟合曲线为： k k k =1 y = Ax 2 + Bx+C, 其系数满足如下正规方程：
N k =1 N k =1 N k =1 N k =1 N k =1 N k =1 N k =1 N k =1 N k =1 N k =1 N N
«数值计算方法»实验报告
③样条插值拟合（三次样条拟合）法结果：
由于三次样条拟合曲线经过每个数据点，故误差 E=0，故只需画出拟合曲线图。 关键代码 1（计算每个子区间上三次样条曲线的系数） ： 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(N)=U(N-1)/B(N-1); for k=N-2:-1:1 M(k+1)=(U(k)-C(k)*M(k+2))/B(k); 实验原理及分析 end M(1)=3*(D(1)-dx0)/H(1)-M(2)/2; 验原理及分析 M(N+1)=3*(dxN-D(N))/H(N)-M(N)/2; 实验原理及分析 for k=0:N-1 S(k+1,1)=(M(k+2)-M(k+1))/6/H(k+1); [X(k+1),X(k+2)]上三次样条曲线 3 次项的系数 S(k+1,2)=M(k+1)/2; %S(k+1,2)为区间 %S(k+1,1)为区间 %M(N+1)的计算公式，具体参看 %M(1)的计算公式，具体参看实 %M(k+1)的求解公式，具体参看
i =0 N N � � 综上所述：N阶Bezier曲线 P(t)=(∑ xi Bi ,N (t), ∑ yi Bi ,N (t )) i =0 i =0
4.实验结果及结论分析：
①线性最小二乘拟合法结果：
A =0.1552,B =59.1848，误差为： E =3.2846 故最小二乘拟合曲线为： y=0.1552x + 59.1848, 拟合误差为： E =3.2846 关键代码 1（求解 A,B 及误差 E） ： xmean=sum(X)/N; tmean=sum(T)/N; B=ymean-A*xmean E=sqrt(sum((A*X+B-T).^2)/N) 关键代码 2（绘图） ： plot(X,A*X+B) hold on plot(X,A*X+B-T,'r') %画出拟合曲线 %在一幅图上画出拟合曲线及误差曲线 %画出误差曲线 %求解时间 X 的平均值 %求解温度 T 的平均值 %系数 B 的计算公式 %计算均方根误差 E
分别利用：线性最小二乘拟合、多项式拟合（最小二乘抛物线拟合 ） 、样条插值 拟合（三次样条拟合） 、三角多项式拟合、贝塞尔曲线拟合这五种方法对这组温 度数据进行拟合，通过拟合的结果比较这五种方法的优缺点。
3.实验原理及分析：
«数值计算方法»实验报告
①线性最小二乘拟合法：
设{( x , y )} 有 N个点，其中横坐标{x } 是确定的。最小二乘拟合曲线为： k k k =1 k k =1 y = Ax + B, 其系数满足如下正规方程：
hk −1 mk −1 + 2(hk −1 + hk ) mk + hk mk −1 = uk , k = 1, 2,..., N − 1
要求出Sk ( x)的最终表达式还差两个端点约束条件，针对三次样条的5种端点约束条件如下：
端点约束条件描述 1.三次紧压样条，确定 S ' ( x0 ), S ' ( x1 )
k
∑x
k =1 N
2 k
− Nx
N
2
∑ (x
k =1 k
− x)
2
其中： x=
∑ xk
k =1
∑y
,y=
k =1
N
N
线性最小二乘法的本质是：多元函数（均方根误差函数）求极值问题。 在处理一些非线性拟合问题时， 我们往往会通过适当的变量代换将问题转化为线 性最小二乘法拟合问题。
②多项式拟合（最小二乘抛物线拟合）法：
«数值计算方法»实验报告
设有N +1个点a{( x j , y j )}N j =0 , 其中y j = f ( x j ), 且横坐标之间等间距，即;
x j = −π +
2 jπ , 其中j = 0,1,..., N N
N k =1
如果f ( x)的周期为2π , 且2M < N , 则存在三角多项式 TM ( x)使得： ∑ ( f ( xk ) − TM ( xk )) 2最小
2.实验内容：
已知某地区一天的温度数据如下： 时间，p.m 1 2 3 4 5 6 7 8 9 10 11 午夜 66 66 65 64 63 63 62 61 60 60 59 58 温度 时间，a.m 1 2 3 4 5 6 7 8 9 10 11 正午 58 58 58 58 57 57 57 58 60 64 67 68 温度
m0 =
包含m0和m N的方程
m 3 m 3 ' (d0 − S' (x0 )) − 1 , mN = (S (xN )− dN−1 )− N−1 h0 2 hN−1 2
2.Natural样条 3.外推样条
m0 = 0, mN = 0 h0 (m2 − m h (m − m ) 1) , mN = mN−1 + N−1 N−1 N−2 h1 hN−2
N
N
其中 : TM ( x) =
⑤贝塞尔曲线拟合法：
贝塞尔曲线具有许多优良的性质，如：任意阶可导、无高阶多项式振荡、 易 于构建等。下面将具体阐述贝塞尔曲线拟合的基本原理：
Bezier曲线可由Bernstein多项式定义，因此有必要提及N 阶Bernstein多项式的定义。 N! i i i N阶Bernstein多项式定义为：Bi , N (t) = CN t (1 − t) N −i , 其中：CN = , i = 0,1, 2,..., N i !( N − i )! N �� � � �� N �� N阶Bezier曲线定义为： P(t) = ∑ PB {P i i , N ( t ), 其中： i }i = 0 为给定的控制点集，即P i =( xi , yi )。
(∑ xk4 ) A + (∑ xk3 ) B + (∑ xk2 )C = ∑ xk2 yk
3 (∑ xk ) A + (∑ xk2 ) B + (∑ xk )C = ∑ xk yk
N k =1