温度分布的曲线拟合

数值计算方法的目的之一是求解一个将自变量与因变量联系起来的拟合函 在科学技术工程和实验中，经常需要从大量的实验数据中寻找拟合曲线，最 标题：温度分布的曲线拟合 比较这五种方法的优缺点（主要考虑误差 «数值计算方法»实验报告 1.1.1.1.实验描述： 简单的是一维情形（一元函数），此时数据的形式为 x 和 y 坐标的有序对，如： 1 1( , ),..., ( , ), { }N Nx y x y x这里的横坐标 是明确的。 数。求解拟合函数的方法有多种，常见的方法有：线性最小二乘拟合、多项式拟 合（最小二乘抛物线拟合）、样条插值拟合（三次样条拟合）、三角多项式拟合、 贝塞尔曲线拟合这五种方法。 本次实验分别利用上述五种方法对一组温度数据进行拟合，通过拟合的结果 ）。 2.2.2.2.实验内容： 已知某地区一天的温度数据如下： 分别利用：线性最小二乘拟合、多项式拟合（最小二乘抛物线拟合）、样条插值 拟合（三次样条拟合）、三角多项式拟合、贝塞尔曲线拟合这五种方法对这组温 度数据进行拟合，通过拟合的结果比较这五种方法的优缺点。 3.3.3.3.实验原理及 分析 定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析 ： 时间，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 1 2 1 1 1 1 1 , ( ) ( ) ( ) {( , )} { } N N k k N N N k k k k k k k N N k k k k N k k k y Ax B x A x B x y x A NB y x y x = = = = = = = = + + = + = ∑ ∑ ∑ ∑ ∑ 设 有 个点，其中横坐标 是确定的。最小二乘拟合曲线为： 其系数满足如下正规方程： 1 1 22 2 1 1 1 1 ( )( ) , ( ) , N N k k k k k k N N k k k k N N k k k k x y N xy x x y y A B y Ax x N x x x x y x y N N = = = = = = − − − = = = − − − = = ∑ ∑ ∑ ∑ ∑ ∑ 解得： 其中： 线性最小二乘法的本质是：多元函数（均方根误差函数）求极值问题。 在处理一些非线性拟合问题时，我们往往会通过适当的变量代换将问题转化为线 性最小二乘法拟合问题。 ②多项式拟合（最小二乘抛物线拟合）法： 1 1 2 4 3 2 2 1 1 1 1 3 2 1 1 1 1 2 1 1 1 + , ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) {( , )} { } N N k k N N N N k k k k k k k k k N N N N k k k k k k k k k N N N k k k k k k N k k k y Ax Bx C x A x B x C x y x A x B x C x y x A x B NC y x y x = = = = = = = = = = = = = = + + + = + + = + + = ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ 设 有 个点，其中横坐标 是确定的。最小二乘抛物线拟合曲线为： 其系数满足如下正规方程： 最小二乘抛物线拟合方法的本质仍然是：多元函数（均方根误差函数）求极 值问题，可见最小二乘法的核心问题是：求解使均方根误差函数取最小值时的拟 合函数的参数。具体方法为（以最小二乘抛物线拟合为例）： 2 2 1 ( , , ) ( ) , , , ( , , ) ( , , ) ( , , ) 0 0, 0 N k k k k E A B C Ax Bx C y A B C E A B C E A B C E A B C A B C = = + + − ∂ ∂ ∂ = = = ∑ 分别对 求偏导数： ， 最终得到上述正规方程。 ③样条插值拟合（三次样条拟合）法： «数值计算方法»实验报告 0 1 2 0 2 3 ,0 ,1 ,2 ,3 1 1 1 ' ' 1 1 1 '' '' 1 1 1 1 ... , ( ) 1. ( ) ( ) ( ) ( ) ( ) 2. ( ) 3. ( ) ( ) 4. ( ) ( ) 5. ( ) ( ) ( ) {( , )} N N k k k k k k k k k k k k k k k k k k k k k k k k N a x x x x b N S x k k S x S x s s x x s x x s x x S x y S x S x S x S x S x S x S x x y = + + + + + + + + + + = < < < < = = = + − + − + − = = = = 设 有 个点，其中: 如果存在 个三次多项式 满足： 则称 '' 1 1 1 3 31 +1 1 0,1,..., 1;2 0,1,..., ;3 5 0,1,..., 2 ( ), , , 0,1,..., 1 6( ), 1, 2,..., 1 ( ) ( ) ( ) ( )( 6 6 6 k k k k k k k k k k k k k k k k k k k k k k k k k N k N k N y y m S x h x x d k N h u d d k N m m y m h S x x x x x x h h h + + − + + = − = − = − − = = − = = − = − = − = − − + − + − 为三次样条函数。（注意：1式中 式中 式中 ） 令 有以上性质得： 1 1 1 1 1 1 ) ( )( ) 6 2( ) , 1, 2,..., 1 ( ) k k k k k k k k k k k k k k y m h x x x h h m h h m h m u k N S x + + − − − − − + − − + + + = = − 要求出 的最终 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 达式还差两个端点约束条件，针对三次样条的5种端点约束条件如下： 端点约束条件描述 0 Nm m包含 和 的方程 ' ' 0 11. ( ), ( )S x S x三次紧压样条，确定 ' ' 11 0 0 0 1 0 1 3 3 ( ( )) , ( ( ) ) 2 2 N N N N N m m m d S x m S x d h h − − − = − − = − − 2.Natural样条 0 0, 0Nm m= = 3.外推样条 0 2 1 1 1 2 0 1 1 1 2 ( ) ( ) , N N N N N N h m m h m m m m m m h h − − − − − − − = − = + 4.抛物线终结样条 0 1 1, N Nm m m m −= = 5.端点曲率调整样条 '' ''0 0( ), ( )N Nm S x m S x= = 其中最常用的是端点约束条件 1：三次紧压样条。 ④三角多项式拟合法： 三角多项式拟合本质上是数字信号处理中的离散余弦变换（DCT），这种方法 可理解为数字信号处理中的 DFT（离散 Fourier 变换）的衍生物，它揭示了事物 内在的周期性。具体的运算方法及公式如下： «数值计算方法»实验报告 0 2 1 1 1 {( , )} , ( ), ; 2 , 0,1,..., ( ) 2 , 2 , ( ) ( ( ) ( )) 2 2 ( )cos( ), ; 0,1,..., ; ( ) sin( ), 1,2,.. N j j j j j j N M k M k k N N j k k j k k k k N a x y y f x j x j N N f x M N T x f x T x a f x jx j M b f x jx j N N π π π = = = = = = − + = < − = = = = ∑ ∑ ∑ 设有 +1个点 其中 且横坐标之间等间距，即 其中 如果 的周期为 且 则存在三角多项式 使得： 最小 其中 其中: .,M 0 1 : ( ) ( cos( ) sin( )) 2 M M j j j a T x a jx b jx = = + +∑其中 ⑤贝塞尔曲线拟合法： 贝塞尔曲线具有许多优良的性质，如：任意阶可导、无高阶多项式振荡、易 于构建等。下面将具体阐述贝塞尔曲线拟合的基本原理： , , 0 0 ! ( ) (1 ) , = , 0,1,2,..., !( )! ( ) ( ), { } =( , ) i i N i i i N N N N N i i N i i i i i i Bezier Bernstein N Bernstein N N Bernstein B t C t t C i N i N i N Bezier P t PB t P P x y − = = = − = − =∑ �� �� �� �� 曲线可由 多项式定义，因此有必要提及 阶 多项式的定义。 阶 多项式定义为： 其中： 阶 曲线定义为： 其中： 为给定的控制点集，即 。 综上所述： , , 0 0 ( )=( ( ), ( )) N N i i N i i N i i N Bezier P t x B t y B t = = ∑ ∑ �� 阶 曲线 4.实验结果及结论分析： ①线性最小二乘拟合法结果： A =0.1552,B =59.1848 E =3.2846 y=0.1552 59.1848, E =3.2846x + ，误差为： 故最小二乘拟合曲线为： 拟合误差为： 关键代码 1（求解 A,B及误差 E）： xmean=sum(X)/N; %求解时间 X 的平均值 tmean=sum(T)/N; %求解温度 T 的平均值 A=sum((X-xmean).*(T-tmean))/sum((X-xmean).^2) %系数 A 的计算公式 B=ymean-A*xmean %系数 B 的计算公式 E=sqrt(sum((A*X+B-T).^2)/N) %计算均方根误差 E 关键代码 2（绘图）： plot(X,A*X+B) %画出拟合曲线 hold on %在一幅图上画出拟合曲线及误差曲线 plot(X,A*X+B-T,'r') %画出误差曲线 «数值计算方法»实验报告 拟合曲线及误差曲线如下： ②多项式拟合（最小二乘抛物线拟合）法结果： 2 A = -0.0564,B =1.5662, C= 53.0707, E =2.2289 y=-0.0564 1.5662 53.0707, E =2.2289x x+ + 误差为： 故最小二乘抛物线拟合为： 拟合误差为： 关键代码 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]; %D 为正规方程的系数矩阵 F=[sum(T.*X.^2);sum(T.*X);sum(T)]; %F 为正规方程的常数向量 P=D\F %求解 A,B,C A=P(1) %A 的值 B=P(2) %B 的值 C=P(3) %C 的值 E=sqrt(sum((A*X.^2+B*X+C-T).^2)/N) %计算均方根误差 E 关键代码 2（绘图）： plot(X,A*X.^2+B*X+C) %画出拟合曲线 hold on %在一幅图上画出拟合曲线及误差曲线 plot(X,A*X.^2+B*X+C-T,'r') %画出误差曲线 拟合曲线及误差曲线如下： «数值计算方法»实验报告 ③样条插值拟合（三次样条拟合）法结果： 由于三次样条拟合曲线经过每个数据点，故误差 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); %M(k+1)的求解公式，具体参看 实验原理及分析 end M(1)=3*(D(1)-dx0)/H(1)-M(2)/2; %M(1)的计算公式，具体参看实 验原理及分析 M(N+1)=3*(dxN-D(N))/H(N)-M(N)/2; %M(N+1)的计算公式，具体参看 实验原理及分析 for k=0:N-1 S(k+1,1)=(M(k+2)-M(k+1))/6/H(k+1); %S(k+1,1)为区间 [X(k+1),X(k+2)]上三次样条曲线 3 次项的系数 S(k+1,2)=M(k+1)/2; %S(k+1,2)为区间 «数值计算方法»实验报告 [X(k+1),X(k+2)]上三次样条曲线 2 次项的系数 S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6; %S(k+1,3)为区间 [X(k+1),X(k+2)]上三次样条曲线 1 次项的系数 S(k+1,4)=T(k+1); %S(k+1,1)为区间 [X(k+1),X(k+2)]上三次样条曲线的常数项 end S ； %S 储存的是每子个区 间上三次样条曲线的系数 关键代码 2（绘图）； plot(X1,T1,X2,T2,X3,T3,X4,T4,X5,T5,X6,T6,X7,T7,X8,T8,X9,T9,X10,T10,X1 1,T11,X12,T12,X13,T13,X14,T14,X15,T15,X16,T16,X17,T17,X18,T18,X19,T19, X20,T20,X21,T21,X22,T22,X23,T23) %plot 命令画出三次样条拟合曲线， 共 23段 三次样条拟合曲线如下： ④三角多项式拟合法结果： [2.1382 1.7128 0.6622 0.0000 - 0.0693 - 0.1230 - 0.1226 - 0.0000 - 0.0081 - 0.0088 0.0740 - 0.0000] [ - 38.5971 - 2.8172 -13.0475 - 3.82 61 - 6.1217 -4.3041 - 2.8762 - A B = = 4.8196 -1.2781 - 5.2767 - 0.5232] 67.6758E = «数值计算方法»实验报告 关键代码 1（计算三角多项式的系数向量 A 和 B 及误差 E）： for j=1:M+1 A(j)=2/N*sum(T0.*cos(j*X0)); %计算 A(j) end for j=1:M B(j)=2/N*sum(T1.*sin(j*X1)); %计算 B(j) end for k=1:N+1 P(k)=A(1)/2+sum(A(1,2:M+1).*cos(X(k)*A0))+sum(B.*sin(X(k)*A0)); %计算 P(k) end A,B %输出 A,B E=sqrt(sum((T-P).^2)/(N+1)) %计算误差 E 关键代码 2（绘图）： plot(X,P,X,P-T,'*r') %画出拟合曲线及误差曲线 拟合曲线及误差曲线如下： «数值计算方法»实验报告 需要再次提及的是：必须先将 X 变换在[-pi,pi]区间上，若缺少这一步则结果 是错误的，这一步非常重要！！！ ⑤贝塞尔曲线拟合法结果： 注意到：Bezier 拟合不宜直接计算误差，故可以从图形上对原始数据曲线图与 拟合曲线图进行比较，从直观上分析其误差。 关键代码： X1=[]; %X1为存储Bezier拟合数据横坐标的向量 T1=X1; %T1为存储Bezier拟合数据纵坐标的向量 for t=0:0.001:1 %t为N阶Bezier曲线的参变量 for i=1:N+1 X0(i)=X(i)*t^(i-1)*(1-t)^(N+1-i)*factorial(N)/factorial(i-1)/factoria l(N+1-i); T0(i)=T(i)*t^(i-1)*(1-t)^(N+1-i)*factorial(N)/factorial(i-1)/factoria l(N+1-i); end X1=[X1 sum(X0)]; %X1为存储Bezier拟合数据横坐标的向量，计 算公式为上面的代码 T1=[T1 sum(T0)]; %T1为存储Bezier拟合数据纵坐标的向量，计 算公式为上面的代码 end plot(X1,T1,X,T,'-*r') %画出Bezier拟合曲线图及原始数曲线据图， 从图形直观上比较误差 Bezier 拟合曲线图及原始数曲线据图如下： «数值计算方法»实验报告 结论分析：通过对上面五种方法的结果进行分析比较，得到如下结论： 1.五种拟合方法的误差由大到小依次为：三角多项式拟合、线性最小二乘拟合、 多项式拟合（最小二乘抛物线拟合）、贝塞尔曲线拟合、样条插值拟合（三次样 条拟合）。 2.线性最小二乘拟合法原理简单，易于求解，但其误差较大的，在拟合精度要求 不高的情况下可优先考虑这种方法，另外许多非线性曲线拟合问题在适当的变换 下可变成线性最小二乘拟合问题，因此：线性最小二乘拟法合适用范围很广，是 一种非常重要的拟合方法。 3.多项式拟合（最小二乘抛物线拟合）法可看成是线性最小二乘拟合法的推广， 实际上两种方法都是通过求解正规方程（具体参见实验原理分析），得到待求系 数，进而得到拟合函数。 4.样条插值拟合（三次样条拟合）法是一种优良的拟合方法，它的拟合误差为零 （因为拟合曲线经过所有的数据点），但这种方法较为复杂，在拟合精度要求较 高时，这种方法无疑是首选。 5.三角多项式拟合法本质上是 DCT（离散余弦变换），在这次实验中其拟合误差 最大，但这并不表示这种方法“一无是处”，事实上，当所研究的问题具有良好 的“周期”性时，这种方法是不二之选。 6.贝塞尔曲线拟合法是非常好的拟合方法，它具有三次样条拟合法的优点，同时 这种方法较为简单、直观，故这种方法应该是较为理想的拟合方法。 附件 1（代码）： ①线性最小二乘拟合法： «数值计算方法»实验报告 X=1:24; %X 为时间坐标向量 T=[58 58 58 58 57 57 57 58 60 64 67 68 66 66 65 64 63 63 62 61 60 60 59 58]; %T 为温度坐标向量 N=24; %N 为数据点数 xmean=sum(X)/N; %求解时间 X 的平均值 tmean=sum(T)/N; %求解温度 T 的平均值 A=sum((X-xmean).*(T-tmean))/sum((X-xmean).^2) %系数 A 的计算公式 B=ymean-A*xmean %系数 B 的计算公式 E=sqrt(sum((A*X+B-T).^2)/N) %计算均方根误差 E plot(X,A*X+B) %画出拟合曲线 hold on %在一幅图上画出拟合曲线及误差曲线 plot(X,A*X+B-T,'r') %画出误差曲线 ②多项式拟合（最小二乘抛物线拟合）法： X=1:24; %X 为时间坐标向量 T=[58 58 58 58 57 57 57 58 60 64 67 68 66 66 65 64 63 63 62 61 60 60 59 58]; %T 为温度坐标向量 N=24; %N 为数据点数 D=[sum(X.^4) sum(X.^3) sum(X.^2);sum(X.^3) sum(X.^2) sum(X);sum(X.^2) sum(X) N]; %D 为正规方程的系数矩阵 F=[sum(T.*X.^2);sum(T.*X);sum(T)]; %F 为正规方程的常数向量 P=D\F %求解 A,B,C A=P(1) %A 的值 B=P(2) %B 的值 C=P(3) %C 的值 E=sqrt(sum((A*X.^2+B*X+C-T).^2)/N) %计算均方根误差 E plot(X,A*X.^2+B*X+C) %画出拟合曲线 hold on %在一幅图上画出拟合曲线及误差曲线 plot(X,A*X.^2+B*X+C-T,'r') %画出误差曲线 ③样条插值拟合（三次样条拟合）法： X=1:24; %X 为时间坐标向量 T=[58 58 58 58 57 57 57 58 60 64 67 68 66 66 65 64 63 63 62 61 60 60 59 58]; %T 为温度坐标向量 N=23; %N 为数据点数减 1，之所以这样取是为了与实验原理中的公式一 致，故此处 N=24-1=23 dx0=-0.8; %dx0=S'(x0),且根据常识凌晨一点钟时温度是随时间增加而减 «数值计算方法»实验报告 小的,故 dx0<0 dxN=-0.5; %dxN=S'(xN),且根据常识午夜零点时温度是随时间增加而减小 的，故 dxN