运用MATLAB解三、四次多项式.doc

MATLAB解多项式
用MATLAB解多项式以三次、四次为例；
运用MATLAB解三、四次多项式，下面以一组数据来演示：
例如，以下数据是某产品从1978年到2010年销售情况，是预测以后几年的销售情况？
1978
1979
1980
1981
1982
1983
1984
1985
1986
566
632
745
755
769
789
985
1110
1313
1987
1988
1989
1990
1991
1992
1993
1994
1995
1428
1782
1920
2150
2292
2601
3149
4338
5145
1996
1997
1998
1999
2000
2001
2002
2003
2004
5809
6241
6854
7656
8772
10007
11374
12567
14332
2005
2006
2007
2008
2009
2010

16614
19228
22844
26404
29688
32074

首先，用MATLAB画图，程序如下：
x=1:33;
y=[566 632 745 755 769 789 985 1110 1313 1428 1782 1920 2150 2292 2601 3149 4338 5145 5809 6241 6854 7656 8772 10007 11374 12567 14332 16614 19228 22844 26404 29688 32074];
plot(x,y)

画出的图形如下图：

近似符合多项式，用多项式进行拟合。
用三次多项式进行拟合：
设三次多项式为y=a0+a1*x+a2*x^2+a3*x^3，用最小二乘法求解多项式系数，程序如下：
format long
x=1:33;
sumx1=sum(x(1,:)); %x中的元素相加%
b=x.*x; %x中各元素的平方%
sumx2=sum(b(1,:)); %x中各元素的平方相加%
c=b.*x; %x中各元素的三次方%
sumx3=sum(c(1,:)); %x中各元素的三次方相加%
d=b.*b; %x中各元素的四次方%
sumx4=sum(d(1,:)); %x中各元素的四次方相加%
e=d.*x; %x中各元素的五次方%
sumx5=sum(e(1,:)); %x中各元素的五次方相加%
f=c.*c; %x中各元素的六次方%
sumx6=sum(f(1,:)); %x中各元素的六次方相加%
A=[33 sumx1 sumx2 sumx3;sumx1 sumx2 sumx3 sumx4;sumx2 sumx3 sumx4 sumx5;sumx3 sumx4 sumx5 sumx6];
y=[566 632 745 755 769 789 985 1110 1313 1428 1782 1920 2150 2292 2601 3149 4338 5145 5809 6241 6854 7656 8772 10007 11374 12567 14332 16614 19228 22844 26404 29688 32074];
sumy1=sum(y(1,:));
i=x.*y;
sumy2=sum(i(1,:));
j=b.*y;
sumy3=sum(j(1,:));
k=c.*y;
sumy4=sum(k(1,:));
Y=[sumy1 sumy2 sumy3 sumy4]';
B=inv(A)*Y

运行结果为：
B =
1.0e+002 *
-5.55549853371922
5.33146148904227
-0.50016139300278
0.01934485988921
即a0=-555.549853371922
a1= 533.146148904227
a2=-50.016139300278
a3=1.934485988921
计算误差，程序如下：
format long
x=1:33;
a=-555.549853371922+533.146148904227.*x-50.016139300278.*x.^2+1.934485988921.*x.^3;
y=[566 632 745 755 769 789 985 1110 1313 1428 1782 1920 2150 2292 2601 3149 4338 5145 5809 6241 6854 7656 8772 10007 11374 12567 14332 16614 19228 22844 26404 29688 32074];
b=a-y;
c=abs(b)./y;
d=sum(c(1,:))/32
结果如下：
d = 0.16981183704177 误差较大。
用四次多项式进行拟合：
设三次多项式为y=a0+a1*x+a2*x^2+a3*x^3，用最小二乘法求解多项式系数，程序如下：
format long
x=1:33;
sumx1=sum(x(1,:)); %x中的元素相加%
b=x.*x; %x中各元素的平方%
sumx2=sum(b(1,:)); %x中各元素的平方相加%
c=b.*x; %x中各元素的三次方%
sumx3=sum(c(1,:)); %x中各元素的三次方相加%
d=b.*b; %x中各元素的四次方%
sumx4=sum(d(1,:)); %x中各元素的四次方相加%
e=d.*x; %x中各元素的五次方%
sumx5=sum(e(1,:)); %x中各元素的五次方相加%
f=c.*c; %x中各元素的六次方%
sumx6=sum(f(1,:)); %x中各元素的六次方相加%
g=c.*d; %x中各元素的七次方%
sumx7=sum(g(1,:)); %x中各元素的七次方相加%
h=d.*d; %x中各元素的八次方%
sumx8=sum(h(1,:)); %x中各元素的八次方相加%
A=[33 sumx1 sumx2 sumx3 sumx4;sumx1 sumx2 sumx3 sumx4 sumx5;sumx2 sumx3 sumx4 sumx5 sumx6;sumx3 sumx4 sumx5 sumx6 sumx7;sumx4 sumx5 sumx6 sumx7 sumx8];
y=[566 632 745 755 769 789 985 1110 1313 1428 1782 1920 2150 2292 2601 3149 4338 5145 5809 6241 6854 7656 8772 10007 11374 12567 14332 16614 19228 22844 26404 29688 32074];
sumy1=sum(y(1,:));
i=x.*y;
sumy2=sum(i(1,:));
j=b.*y;
sumy3=sum(j(1,:));
k=c.*y;
sumy4=sum(k(1,:));
l=d.*y;
sumy5=sum(l(1,:))
Y=[sumy1 sumy2 sumy3 sumy4 sumy5]';
B=inv(A)*Y

运行结果：
B =
1.0e+003 *
1.05813149711536
-0.30415440597106
0.05699047661168
-0.00291140395249
0.00007126308737
即a0=1058.13149711536
a1=-304.15440597106
a2=56.99047661168
a3=-2.91140395249
a4=0.07126308737
计算误差，程序如下：
format long
x=1:33;
a=1058.13149711536-304.15440597106.*x+56.99047661168*x.^2-2.91140395249.*x.^3+.07126308737.*x.^4;
y=[566 632 745 755 769 789 985 1110 1313 1428 1782 1920 2150 2292 2601 3149 4338 5145 5809 6241 6854 7656 8772 10007 11374 12567 14332 16614 19228 22844 26404 29688 32074];
b=a-y;
c=abs(b)./y;
d=sum(c(1,:))/32
结果如下：
d = 0.08262969345569
用图形观察拟合效果：
程序如下：
ezplot('1058.13149711536-304.15440597106*x+56.99047661168*x^2-2.91140395249*x^3+.07126308737*x^4',1:33);
hold on
x=1:33;
y=[566 632 745 755 769 789 985 1110 1313 1428 1782 1920 2150 2292 2601 3149 4338 5145 5809 6241 6854 7656 8772 10007 11374 12567 14332 16614 19228 22844 26404 29688 32074];
plot(x,y,'r-o')
图形如下：

可看出拟合效果很好。
用四次进行拟合预测
方程为：
y=058.13149711536-304.15440597106.*x+56.99047661168*x.^2-2.91140395249.*x.^3+.07126308737.*x.^4
将2011年化成34带入预测，即预测2011年经x=34带入既得结果，以后年份以此类推。