数值分析上机实验报告选题：曲线拟合的最小二乘法指导老师：专业：学号：姓名：课题八 曲线拟合的最小二乘法一、问题提出从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。

在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量y 与时间t 的拟合曲线。

二、要求1、用最小二乘法进行曲线拟合；2、近似解析表达式为()33221t a t a t a t ++=ϕ；3、打印出拟合函数()t ϕ，并打印出()j t ϕ与()j t y 的误差，12,,2,1 =j ；4、另外选取一个近似表达式，尝试拟合效果的比较；5、*绘制出曲线拟合图*。

三、目的和意义1、掌握曲线拟合的最小二乘法；2、最小二乘法亦可用于解超定线代数方程组；3、探索拟合函数的选择与拟合精度间的关系。

四、计算公式对于给定的测量数据(x i ,f i )(i=1,2,…，n )，设函数分布为∑==mj j j x a x y 0)()(ϕ特别的，取)(x j ϕ为多项式j j x x =)(ϕ (j=0, 1,…，m )则根据最小二乘法原理，可以构造泛函∑∑==-=ni mj i j j i m x a f a a a H 110))((),,,(ϕ令0=∂∂ka H(k=0, 1,…，m ) 则可以得到法方程⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100m m m m m m m m f f f a a a ϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕ求该解方程组，则可以得到解m a a a ,,,10 ，因此可得到数据的最小二乘解∑=≈mj j j x a x f 0)()(ϕ曲线拟合：实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。

曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。

五、结构程序设计在程序结构方面主要是按照顺序结构进行设计，在进行曲线的拟合时，为了进行比较，在程序设计中，直接调用了最小二乘法的拟合函数polyfit ，并且依次调用了plot 、figure 、hold on 函数进行图象的绘制，最后调用了一个绝对值函数abs 用于计算拟合函数与原有数据的误差，进行拟合效果的比较。

5.1用一元三次多项式()33221t a t a t a t ++=ϕ进行拟合计算解析表达式系数： a1, a2, a3t=[0 5 10 15 20 25 30 35 40 45 50 55];y=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64]; >> n=length(xi);f=0.34364.*10.^(-4)*x.^3-5.2156.*10.^(-3)*x.^2+0.26340.*x+0.017839;x=0:0.01:55;F=0.34364.*10.^(-4)*x.^3-5.2156.*10.^(-3)*x.^2+0.26340.*x+0.017839;fy=abs(f-y);fy2=fy.^2;Ew=max(fy),E1=sum(fy)/n,E2=sqrt((sum(fy2))/n)plot(xi,y,'t*'), hold on, plot(t,F,'b-'), hold off所得函数为4332(t)0.3436410 5.2156100.26340.013839t t t ϕ--=⨯-⨯++ 运行后屏幕显示数据),(i i y x 与拟合函数f 的最大误差Ew ，平均误差E1和均方根误差E2及其数据点),(i i y x 和拟合曲线y=f(x)的图形如图5.1.Ew =0.4243E1 =0.0911E2 =0.1467图5.1一元三次多项式拟合曲线误差图5.2用一元四次多项式()4433221t a t a t a t a t +++=ϕ进行拟合：计算多项式系数：a 1, a 2, a 3, a 4xi=[0 5 10 15 20 25 30 35 40 45 50 55];y=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64];n=length(xi);x=0:0.01:55;f=0.6026.*10.^(-6)*x.^4-0.31918.*10.^(-4)*x.^3-0.0029323.*x.^2+0.23807.*x+0.060449;x=0:0.01:55;F=0.6026.*10.^(-6)*x.^4-0.31918.*10.^(-4)*x.^3-0.0029323.*x.^2+0.23807.*x+0.060449;fy=abs(f-y);fy2=fy.^2;Ew=max(fy),E1=sum(fy)/n,E2=sqrt((sum(fy2))/n) plot(xi,y,'r*'), hold on, plot(x,F,'b-'), hold off所得函数为644332(t)0.6026100.3191810 2.9323100.238070.060449t t t t ϕ---=⨯-⨯-⨯++运行后屏幕显示数据),(i i y x 与拟合函数f 的最大误差Ew ，平均误差E1和均方根误差E2及其数据点),(i i y x 和拟合曲线y=f(x)的图形如图5.2。

Ew = 0.3897 E1 = 0.1034、 E2 =0.1429图5.2一元四次多项式拟合曲线误差图5.3用一元二次多项式()2321t a t a a t ++=ϕ进行拟合：计算多项式系数：a 1, a 2, a 3输入程序： >> syms a1 a2 a3x=[0 5 10 15 20 25 30 35 40 45 50 55]; fi=a1.*x.^2+ a2.*x+ a3运行后屏幕显示关于a 1,a 2和a 3的线性方程组：fi=[ a3, 25*a1 + 5*a2 + a3, 100*a1 + 10*a2 + a3, 225*a1 + 15*a2 + a3, 400*a1 + 20*a2 + a3, 625*a1 + 25*a2 + a3, 900*a1 + 30*a2 + a3, 1225*a1 + 35*a2 + a3, 1600*a1 + 40*a2 + a3, 2025*a1 + 45*a2 + a3, 2500*a1 + 50*a2 + a3, 3025*a1 + 55*a2 + a3]编写构造误差平方和的MATLAB 程序：y=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64];fi =[ a3, 25*a1 + 5*a2 + a3, 100*a1 + 10*a2 + a3, 225*a1 + 15*a2 + a3, 400*a1 + 20*a2 + a3, 625*a1 + 25*a2 + a3, 900*a1 + 30*a2 + a3, 1225*a1 + 35*a2 + a3, 1600*a1 + 40*a2 + a3, 2025*a1 + 45*a2 + a3, 2500*a1 + 50*a2 + a3, 3025*a1 + 55*a2 + a3];fy=fi-y;fy2=fy.^2;J=sum(fy.^2) 运行后屏幕显示误差平方和如下：J =(100*a1 + 10*a2 + a3 - 54/25)^2 + (25*a1 + 5*a2 + a3 - 127/100)^2 + (225*a1 + 15*a2 + a3 - 143/50)^2 + (400*a1 + 20*a2 + a3 - 86/25)^2 + (900*a1 + 30*a2 + a3 - 83/20)^2 + (625*a1 + 25*a2 + a3 - 387/100)^2 + (1225*a1 + 35*a2 + a3 - 437/100)^2 + (1600*a1 + 40*a2 + a3 - 451/100)^2 + (2025*a1 + 45*a2 + a3 - 229/50)^2 + (2500*a1 + 50*a2 + a3 - 201/50)^2 + (3025*a1 + 55*a2 + a3 - 116/25)^2 + a3^2为求4321,,,a a a a 使J 达到最小，只需利用极值的必要条件0=∂∂ka J)4,3,2,1(=k ，得到关于4321,,,a a a a 的线性方程组，这可以由下面的MATLAB 程序完成，即输入程序 ：>> syms a1 a2 a3J =(100*a1 + 10*a2 + a3 - 54/25)^2 + (25*a1 + 5*a2 + a3 - 127/100)^2 + (225*a1 + 15*a2 + a3 - 143/50)^2 + (400*a1 + 20*a2 + a3 - 86/25)^2 + (900*a1 + 30*a2 + a3 - 83/20)^2 + (625*a1 + 25*a2 + a3 - 387/100)^2 + (1225*a1 + 35*a2 + a3 - 437/100)^2 + (1600*a1 + 40*a2 + a3 - 451/100)^2 + (2025*a1 + 45*a2 + a3 - 229/50)^2 + (2500*a1 + 50*a2 + a3 - 201/50)^2 + (3025*a1 + 55*a2 + a3 - 116/25)^2 + a3^2;Ja1=diff(J,a1);Ja2=diff(J,a2);Ja3=diff(J,a3);Ja11=simple(Ja1),Ja21=simple(Ja2),Ja31=simple(Ja3), 运行后屏幕显示J 分别对a1, a2 ,a3的偏导数如下： Ja11 =49967500*a1 + 1089000*a2 + 25300*a3 - 217403/2 Ja21 = 1089000*a1 + 25300*a2 + 660*a3 - 27131/10 Ja31 = 25300*a1 + 660*a2 + 24*a3 - 3987/50解线性方程组Ja11 =0，Ja21 =0，Ja31 =0输入下列程序： >>A=[49967500,1089000,25300;1089000, 25300,660;25300,660,24]; B=[217403/2,27131/10,3987/50]; C=B/A, F=poly2sym(C)运行后屏幕显示拟合函数f 及其系数C 如下： C =-0.0024 0.2037 0.2305F = (7338734818964133*x)/36028797018963968 - (5489104202452799*x^2)/2305843009213693952 + 8303449950332545/36028797018963968故所求的拟合曲线为：0.230470.203690.0023805)(22++-=x x x f编写下面的MATLAB 程序估计其误差，并作出拟合曲线和数据的图形。