77翟又文上机作业.doc

昆明理工大学 《数值分析》上机实验报告 姓名 ： 翟又文 学号 ： 2013230077 专业 ： 材料加工工程 学院 ： 材料科学与工程学院 授课教师 ： 李玉兰 时间 ： 2013.12.20 一、课题名称：曲线拟合的最小二乘法 二、问题的提出 对测量数据的拟合是一个既古老，但又非常实用的问题。所谓的拟合就是插值与逼近的统称。给定一组杂乱无章的实验数据，其中，要求构造一条曲线顺序通过这些数据点，称为对这些数据点进行插值(interpolation),所构造的曲线称为插值曲线。在某些情况下，实验所得到的数据点本身就比较粗糙，要求构造一条曲线严格通过给定的一组数据点没有什么实际意义，更合理的方法是构造一条曲线使之在某种意义下最为接近给定的数据点，这种方法称之为对已知数据点的逼近(approximation),所构造的曲线称为逼近曲线。 对于具体的实验，通常不是先给出函数的解析式，然后再进行实验，而是通过对实验的观察和测量给出一定数量的离散的数据点，再根据得到的这些数据点来求出具体的函数解析式。又因为实验中存在着诸如测量误差等的误差，这就导致实际真实的解析式曲线并不一定通过测量得到的所有数据点。最小二乘法是求解这一问题很好的方法，故本实验采用最小二乘法来实现对所给数据点的拟合。 三、实验目的和意义 1、掌握曲线拟合的最小二乘法； 2、最小二乘法亦可用于解超定线代数方程组； 3、探索拟合函数的选择与拟合精度间的关系。 四、数学原理 最小二乘法的一般提法是对于给定的测量数据，要求在给定函数类中找到一个函数，使满足： 比较一般的要求是需满足下式： 其中的表达式如下： 也就是说，要对已知的数据点进行拟合，就需要求出函数，而求函数的关键就是求出函数的系数，使得下式取得极小值。 引入内积记号 则可得到如下式所示的法方程。 其中在上线性无关。 设的任意线性组合在点集上至多只有个不同的零点，则称在点集上满足条件。 很显然在任意个点上满足条件。 因此，如果在上满足条件，则有 求出函数的系数，就能得到函数的表达式，并可以证明 由此，可知函数是最小二乘解。 实际应用中常常使用多项式拟合。即取 此时拟合函数为： 法方程为 对于其他形式的拟合函数，通常将其转换为线性模型，即转换为如下式所示的模型。 例如，对于形如的拟合函数，我们通常先对其取对数可得 然后令，，，得到如下所示的线性模型 求和，且满足 解法方程组得到和后，由下式求得和。 求出和后将其代入拟合函数可得 五、实验内容及要求 1、实验内容 从随机的数据中找出其规律性，给处近似表达式的问题，在生产和实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中，根据统计数据的含碳量与时间的关系，试求含碳量和时间的拟合曲线。 含碳量和时间的关系如下表所示。 (分) 0 5 10 15 20 25 30 35 40 45 50 55 () 0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64 2、实验要求 （1）用最小二乘法进行曲线拟合； （2）近似解析表达式为； （3）打印出拟合函数，并打印出与的误差，； （4）另外选取一个近似表达式，尝试拟合效果的比较； （5）绘制出曲线拟合图。 六、计算过程 方法1 解：设近似解析表达式为，则有，，，，，即，计算离散内积有 解法方程组 得，， 因此近似解析表达式为 各点的拟合误差为： 拟合平方误差为： 拟合均方误差为： 拟合误差的最大值为： 曲线拟合图如图(1)所示 图 (1) 采用近似解析表达式 拟合时，各数据点的拟合误差以及拟合平方误差和拟合均方误差如表(1)所示： 表(1) 序号 t y 拟合值 拟合误差 平方误差 均方误差 1 0 0 0 0 1.13 1.06 2 5 1.27 1.18 7493 3 10 2.16 2.14 514 4 15 2.86 2.89 758 5 20 3.44 3.46 398 6 25 3.87 3.88 108 7 30 4.15 4.17 615 8 35 4.37 4.37 1.03 9 40 4.51 4.49 448 10 45 4.58 4.56 390 11 50 4.62 4.61 120 12 55 4.64 4.66 447 由表(1)可知采用方法1得到的拟合曲线其拟合误差在处最大。 方法2 解：设近似解析表达式为其中和为正常数，对该解析表达式等号两边取对数可得，令，则有 所以，计算离散内积有 解法方程组 得，， 综上所述，近似解析表达式为 各点的拟合误差为： 拟合平方误差为： 拟合均方误差为： 拟合误差的最大值为： 曲线拟合图如图(2)所示。 图 (2) 采用近似解析表达式 进行拟合时，各数据点的拟合误差以及拟合平方误差和拟合均方误差如表(2)所示： 表 (2) 序号 t y 拟合值 拟合误差 平方误差 均方误差 1 0 0 0 0 3.37 5.81 2 5 1.27 1.16 1112 3 10 2.16 2.46 9265 4 15 2.86 3.16 9237 5 20 3.44 3.58 2101 6 25 3.87 3.86 3.556 7 30 4.15 4.06 773 8 35 4.37 4.21 2570 9 40 4.51 4.32 3464 10 45 4.58 4.41 2727 11 50 4.62 4.49 1715 12 55 4.64 4.55 798 由表(1)可知采用方法1得到的拟合曲线其拟合误差在处最大。 七、结构程序设计 方法1源程序为： clear %清除变量 clc %清除命令窗口 clf %清除图形窗口 %给定测量数据点(t,y) t=[0 5 10 15 20 25 30 35 40 45 50 55] y=[0 1.27*10^(-4) 2.16*10^(-4) 2.86*10^(-4) 3.44*10^(-4) 3.87*10^(-4) 4.15*10^(-4) 4.37*10^(-4) 4.51*10^(-4) 4.58*10^(-4) 4.62*10^(-4) 4.64*10^(-4)] %计算给定的数据点数目n n=length(y) %给定需要拟合的数据的最高次多项式的次数 m=3 for k=1:m g=zeros(1,m) %构造1行三列矩阵 for j=1:m nt=0 for i=1:n nt=nt+(t(i).^j).*(t(i).^k) %求内积 end g(j)=nt end A(k,:)=g %求系数矩阵 ny=0 for i=1:n ny=ny+y(i).*(t(i).^k) end b(k,1)=ny end a=A\b t1=[0:5:55] f=a(1)*t1+a(2)*t1.^2+a(3)*t1.^3 %求拟合曲线 plot(t1,f,':b') %画出拟合曲线 grid on %画网格线 hold on %保留拟合曲线 plot(t(:),y(:),'*r') %画出测量数据点 hold on %保留测量数据点 plot(t,y,'r') %画测量数据点曲线 hold on legend('拟合曲线','测量数据点','测量数据点曲线',2) %添加图例 %求拟合平方误差 w1=0 for i=1:12 w(i)=(f(i)-y(i))^2 w1=w1+w(i) end %求均方误差 jw1=w1.^0.5 %输出函数y的表达式、拟合平方误差及拟合均方误差 fprintf('方法1解析表达式为：y=%11.10f*t+%11.10f*t^2+%11.10f*t^3\n',a(1),a(2),a(3)) fprintf('方法1的拟合平方误差为：%13.13f\n',w1) fprintf('方法1的拟合均方误差为：%13.13f',jw1) 方法2源程序为： clear %清除变量 clc %清除命令窗口 clf %清除图形窗口 %给定测量数据点(t,y) t=[0 5 10 15 20 25 30 35 40 45 50 55] y=[0 1.27*10^(-4) 2.16*10^(-4) 2.86*10^(-4) 3.44*10^(-4) 3.87*10^(-4) 4.15*10^(-4) 4.37*10^(-4) 4.51*10^(-4) 4.58*10^(-4) 4.62*10^(-4) 4.64*10^(-4)] %计算给定的数据点数目n n=length(y) %给定系数矩阵的行数 m=2 for k=1:m g=zeros(1,m) %构造1行2列矩阵 for j=1:m nt=0 for i=2:n if k==1&j==1 nt=nt+1 elseif k==1&j==2|k==2&j==1 nt=nt-1/t(i) elseif k==2&j==2 nt=nt+1/(t(i).^2) end end g(j)=nt end A(k,:)=g ny=0 for i=2:n if k==1&j==2 ny=ny+log(y(i)) elseif k==2&j==2 ny=ny-log(y(i))*(1/t(i)) end end b(k,1)=ny end a1=A\b a=exp(a1(1)) t1=[0:5:55] f=a*exp(-a1(2)./t1) plot(t1,f,'--b') grid on %画网格线 hold on %保留拟合曲线 plot(t(:),y(:),'*r') %画出测量数据点 hold on plot(t,y,'r') hold on legend('拟合曲线','测量数据点','测量数据点曲线',2) %添加图例 %求拟合平方差 w1=0 for i=1:12 w(i)=(f(i)-y(i))^2 w1=w1+w(i) end %求均方误差 jw1=w1.^0.5 %输出函数y的表达式、拟合平方误差及拟合均方误差 fprintf('方法2的解析表达式为：y=%f*e^(-%f/t)\n',a,a1(2)) fprintf('方法2的拟合平方误差为：%13.13f\n',w1) fprintf('方法2的拟合均方误差为：%13.13f',jw1) 拟合结果比较 clear %清除变量 clc %清除命令窗口 clf %清除图形窗口 %给定测量数据点(t,y) t=[0 5 10 15 20 25 30 35 40 45 50 55] y=[0 1.27*10^(-4) 2.16*10^(-4) 2.86*10^(-4) 3.44*10^(-4) 3.87*10^(-4) 4.15*10^(-4) 4.37*10^(-4) 4.51*10^(-4) 4.58*10^(-4) 4.62*10^(-4) 4.64*10^(-4)] %计算给定的数据点数目n n=length(y) %给定需要拟合的数据的最高次多项式的次数 m=3 for k=1:m g=zeros(1,m) %构造1行三列矩阵 for j=1:m nt=0 for i=1:n nt=nt+(t(i).^j).*(t(i).^k) %求内积 end g(j)=nt end A(k,:)=g %求系数矩阵 ny=0 for i=1:n ny=ny+y(i).*(t(i).^k) end b(k,1)=ny end a=A\b t1=[0:5:55] f=a(1)*t1+a(2)*t1.^2+a(3)*t1.^3 %求拟合曲线 plot(t1,f,':b') %画出拟合曲线 grid on %画网格线 hold on %保留拟合曲线 plot(t(:),y(:),'*r') %画出测量数据点 hold on %保留测量数据点 plot(t,y,'r') hold on %求拟合平方误差 w1=0 for i=1:12 w(i)=(f(i)-y(i))^2 w1=w1+w(i) end %求均方误差 jw1=w1.^0.5 %输出函数y的表达式、拟合平方误差及拟合均方误差 fprintf('方法1解析表达式为：y=%11.10f*t+%11.10f*t^2+%11.10f*t^3\n',a(1),a(2),a(3)) fprintf('方法1的拟合平方误差为：%13.13f\n',w1) fprintf('方法1的拟合均方误差为：%13.13f',jw1) clear %清除变量 clc %清除命令窗口 %给定测量数据点(t,y) t=[0 5 10 15 20 25 30 35 40 45 50 55] y=[0 1.27*10^(-4) 2.16*10^(-4) 2.86*10^(-4) 3.44*10^(-4) 3.87*10^(-4) 4.15*10^(-4) 4.37*10^(-4) 4.51*10^(-4) 4.58*10^(-4) 4.62*10^(-4) 4.64*10^(-4)] %计算给定的数据点数目n n=length(y) %给定系数矩阵的行数 m=2 for k=1:m g=zeros(1,m) %构造1行2列矩阵 for j=1:m nt=0 for i=2:n if k==1&j==1 nt=nt+1 elseif k==1&j==2|k==2&j==1 nt=nt-1/t(i) elseif k==2&j==2 nt=nt+1/(t(i).^2) end end g(j)=nt end A(k,:)=g ny=0 for i=2:n if k==1&j==2 ny=ny+log(y(i)) elseif k==2&j==2 ny=ny-log(y(i))*(1/t(i)) end end b(k,1)=ny end a1=A\b a=exp(a1(1)) t1=[0:5:55] f=a*exp(-a1(2)./t1) plot(t1,f,'--b') grid on %画网格线 hold on %保留拟合曲线 plot(t(:),y(:),'*r') %画出测量数据点 hold on plot(t,y,'r') hold on legend('法1拟合曲线','测量数据点','测量数据点曲线','法2拟合曲线',2) %添加图例 %求拟合平方差 w1=0 for i=1:12 w(i)=(f(i)-y(i))^2 w1=w1+w(i) end %求均方误差 jw1=w1.^0.5 %输出函数y的表达式、拟合平方误差及拟合均方误差 fprintf('方法2的解析表达式为：y=%f*e^(-%f/t)\n',a,a1(2)) fprintf('方法2的拟合平方误差为：%13.13f\n',w1) fprintf('方法2的拟合均方误差为：%13.13f',jw1) 八、结果分析和讨论 方法1运行结果截图如图(3)所示。 图 (3) 方法2运行结果截图如图(4)所示。 图(4) 两种拟合方法得到的曲线比较结果如图(5)所示。 图 (5) 由表（1）和表（2）可知，采用方法1所述的近似解析表达式进行曲线拟合得到的均方误差为，而选用方法二所述的近似解析表达式进行曲线的拟合时得到的均方误差为，误差近似为误差的5.4687倍，再结合图(5)，可以清楚的看出采用方法1所述的近似解析表达式进行曲线拟合时拟合精度比采用方法2所述的近似解析表达式好，对于相同的实验数据采用的表达式不同，拟合精度也不相同。