曲线拟合及其MATLAB 程序
实验任务：用最小二乘法将一组离散数据进行曲线拟合，并用MATLAB 做出曲线。
实验要求：
第三次作业的最后一题求函数关系并作出曲线，此组x 和y 的实验数据如下表所示。 x i
1 2 3 4 5 6 7 8 9 10 y i 6.00 3.50 2.67 2.75 2.00 1.80 1.70 1.63 1.56 1.50
实验内容：
（1）在在MATLAB 工作窗口输入程序
>> x=[1 2 3 4 5 6 7 8 9 10];
y=[6.00 3.50 2.67 2.75 2.00 1.80 1.70 1.63 1.56 1.50];
plot(x,y,'r*'),
legend('实验数据(xi,yi)')
xlabel('x'), ylabel('y'),
title('实验数据点(xi,yi)的散点图')
（2）运行后屏幕显示数据的散点图（略）。
（3）编写下列MATLAB 程序计算)(x f 在),(i i y x 处的函数值，即输入程序 >> syms a1 a2 a3 a4;
x=[1 2 3 4 5 6 7 8 9 10];
fi=a1.*x.^3+ a2.*x.^2+ a3.*x+ a4
运行后屏幕显示关于a1, a2, a3和a4的线性方程组
fi =
[a1+a2+a3+a4, 8*a1+4*a2+2*a3+a4, 27*a1+9*a2+3*a3+a4, 64*a1+16*a2+4*a3+a4, 125*a1+25*a2+5*a3+a4,
216*a1+36*a2+6*a3+a4, 343*a1+49*a2+7*a3+a4,
512*a1+64*a2+8*a3+a4, 729*a1+81*a2+9*a3+a4,
1000*a1+100*a2+10*a3+a4]
编写构造误差平方和的MATLAB 程序
>> y=[6.00 3.50 2.67 2.75 2.00 1.80 1.70 1.63 1.56 1.50];
fi
=[a1+a2+a3+a4,8*a1+4*a2+2*a3+a4,27*a1+9*a2+3*a3+a4,64*a1+16*a2+4*a3+a4,125*a1+25*a2+5*a3+a4,216*a1+36*a2+6*a3+a4,343*a1+49*a2+7*a3+a4,512*a1+64*a2+8*a3+a4,729*a1+81*a2+9*a3+a4, 1000*a1+100*a2+10*a3+a4]; fy=fi-y; fy2=fy.^2;J=sum(fy.^2)
运行后屏幕显示误差平方和如下
J =
(a1+a2+a3+a4-6)^2+(8*a1+4*a2+2*a3+a4-7/2)^2+(27*a1+9*a2+3*a3+a4-267/100)^2+(64*a1+16*a2+4*a3+a4-11/4)^2+(125*a1+25*a2+5*a3+a4-2)^2+(216*a 1+36*a2+6*a3+a4-9/5)^2+(343*a1+49*a2+7*a3+a4-17/10)^2+(512*a1+64*a2+8*a3+a4-163/100)^2+(729*a1+81*a2+9*a3+a4-39/25)^2+(1000*a1+100*a2+10*a3+a4-3/2)^2
为求4321,,,a a a a 使J 达到最小，只需利用极值的必要条件0
=??k a J 。)4,3,2,1(=k ，得到关于4321,,,a a a a 的线性方程组，这可以由下面的MATLAB 程序完成，即输入程序
>> syms a1 a2 a3 a4
J=(a1+a2+a3+a4-6)^2+(8*a1+4*a2+2*a3+a4-7/2)^2+(27*a1+9*a2+3*a3+a4
-267/100)^2+(64*a1+16*a2+4*a3+a4-11/4)^2+(125*a1+25*a2+5*a3+a4-2)^2+( 216*a1+36*a2+6*a3+a4-9/5)^2+(343*a1+49*a2+7*a3+a4-17/10)^2+(512*a1+6 4*a2+8*a3+a4-163/100)^2+(729*a1+81*a2+9*a3+a4-39/25)^2+(1000*a1+100* a2+10*a3+a4-3/2)^2;
Ja1=diff(J,a1); Ja2=diff(J,a2); Ja3=diff(J,a3); Ja4=diff(J,a4);
Ja11=simple(Ja1), Ja21=simple(Ja2), Ja31=simple(Ja3), Ja41=simple(Ja4),
运行后屏幕显示J分别对a1, a2 ,a3 ,a4的偏导数如下
Ja11 =6050*a4+3956810*a1+441650*a2+50666*a3-497579/50
Ja21 =770*a4+441650*a1+50666*a2+6050*a3-66681/50
Ja31 =110*a4+50666*a1+6050*a2+770*a3-10679/50
Ja41 = 20*a4+6050*a1+770*a2+110*a3-2511/50
解线性方程组Ja11 =0，Ja21 =0，Ja31 =0，Ja41 =0，输入下列程序
>>A=[3956810, 441650, 50666, 6050; 441650, 50666, 6050, 770; 50666, 6050, 770, 110; 6050, 770, 110, 20];
B=[497579/50, 66681/50, 10679/50, 2511/50];
C=B/A, f=poly2sym(C)
运行后屏幕显示拟合函数f及其系数C如下
C = -0.0176 0.3744 -2.6452 7.9627
f=-5066427105276917/288230376151711744*x^3+6744376644764349/1801439 8509481984*x^2-2978180218619021/1125899906842624*x+112064********* 3/140737488355328
故所求的拟合曲线为
f(x)= -0.0176 x3 + 0.3744 x2 -2.6452 x+ 7.9627 .
（4）编写下面的MATLAB程序估计其误差，并作出拟合曲线和数据的图形，输入程序
>> xi =[1 2 3 4 5 6 7 8 9 10];
y=[6.00 3.50 2.67 2.75 2.00 1.80 1.70 1.63 1.56 1.50];
n=length(xi);
f= -0.0176.*xi.^3+0.3744.*xi.^2-2.6452.*xi +7.9627;
x=1: 1: 10;
F= -0.0176.*xi.^3+0.3744.*xi.^2 -2.6452.*xi +7.9627; 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 legend('数据点(xi,yi)','拟合曲线y=f(x)'), xlabel('x'), ylabel('y'),
title('实验数据点(xi,yi)和拟合曲线y=f(x)的图形') 运行后屏幕显示数据),(i i y x 与拟合函数f 的最大误差Ew ，平均误差E1和均方根误差E2及其数据点),(i i y x 和拟合曲线y=f(x)的图形如下所示。 Ew = 0.5291
E1 = 0.2163
E2 = 0.2763
