曲线拟合的最小二乘法matlab举例

曲线拟合的最小二乘法

学院：光电信息学院 姓名：赵海峰

学号：200820501001

一、曲线拟合的最小二乘法原理：

由已知的离散数据点选择与实验点误差最小的曲线

)(...)()()(1100x a x a x a x S n n ϕϕϕ+++=

称为曲线拟合的最小二乘法。 若记

),()()(),(0

i k i j m

i i k j x x x ϕϕωϕϕ∑==

k i k i m

i i k d x x f x f ≡=∑=)()()(),(0

ϕωϕ

上式可改写为),...,1,0(;),(n k d a k j n

o

j j k -=∑=ϕϕ这个方程成为法方程，可写成距阵

形式

d Ga =

其中,),...,,(,),...,,(1010T n T n d d d d a a a a ==

⎥⎥⎥⎥

⎦⎤⎢⎢⎢⎢⎣⎡=),(),(),()(),(),(),(),(),(10

1110101000n n n n n n G ϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕϕ 。

它的平方误差为：.)]()([)(||||20

22i i m

i i x f x S x -=

∑=ωδ

二、数值实例：

下面给定的是乌鲁木齐最近1个月早晨7:00左右(新疆时间)的天气预报所得

到的温度数据表，按照数据找出任意次曲线拟合方程和它的图像。

下面应用Matlab编程对上述数据进行最小二乘拟合

三、Matlab程序代码：

x=[1:1:30];

y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1]; a1=polyfit(x,y,3) %三次多项式拟合%

a2= polyfit(x,y,9) %九次多项式拟合%

a3= polyfit(x,y,15) %十五次多项式拟合%

b1=polyval(a1,x)

b2=polyval(a2,x)

b3=polyval(a3,x)

r1= sum((y-b1).^2) %三次多项式误差平方和%

r2= sum((y-b2).^2) %九次次多项式误差平方和%

r3= sum((y-b3).^2) %十五次多项式误差平方和%

plot(x,y,'*') %用*画出x,y图像%

hold on

plot(x,b1, 'r') %用红色线画出x,b1图像%

hold on

plot(x,b2, 'g') %用绿色线画出x,b2图像%

hold on

plot(x,b3, 'b:o') %用蓝色o线画出x,b3图像%

四、数值结果：

不同次数多项式拟和误差平方和为：

r1 = 67.6659

r2 = 20.1060

r3 = 3.7952

r1、r2、r3分别表示三次、九次、十五次多项式误差平方和。

拟和曲线如下图：

上图中*代表原始数据，红色曲线代表三次多项式拟合曲线，绿色曲线代表九次多项式拟合曲线，蓝色o线代表十五次多项式拟合曲线。

五、结论：

以上结果可以看到用最小二乘拟合来求解问题时，有时候他的结果很接近实际情况，有时候跟实际情况里的太远，因为所求得多项式次数太小时数据点之间差别很大，次数最大是误差最小但是有时后不符合实际情况，所以用最小二乘法时次数要取合适一点。

从上面的拟合中也可以得到多项式拟合误差平方和随着拟合多项式次数的增加而逐渐减小，拟合的曲线更靠近实际数据。拟合更准确。

如：根据数据表，用最小二乘法求三次多项式拟合这组数据。数据表

x 0

1 2 3 4 5 6 7 8 y 1 2.5 2 1 3 5 5.5 4

3

#include "stdio.h" #include "conio.h" #include "stdlib.h" #include "math.h" #define N 9 #define M 2 #define K 2*M

void zhuyuan (int k,int n,float a[M+1][M+2]) { int t,i,j; float x,y;

x=fabs(a[k][k]);t=k; for (i=k+1;i<=n;i++) if (fabs(a[i][k])>x) {x=fabs(a[i][k]);t=i;} for (j=k;j<= n+1;j++)

{y=a[k][j];a[k][j]=a[t][j];a[t][j]=y;} }

void xiaoyuan(int n,float a[M+1][M+2]) { int k,i,j;

for(i=0;i

for (j=i+1;j<=n;j++) for (k=i+1;k<=n+1;k++)

a[j][k]=a[j][k]-a[j][i]*a[i][k]/a[i][i]; } }

void huidai(int n,float a[M+1][M+2],float x[M+1]) { int i,j;

x[n]=a[n][n+1]/a[n][n]; for (i=n-1;i>=0;i--) { x[i]=a[i][n+1]; for (j=i+1;j<=n;j++) x[i]=x[i]-a[i ][j]*x[j]; x[i]=x[i]/a[i][i]; } }

wk_ad_begin({pid : 21});wk_ad_after(21, function(){$('.ad-hidden').hide();}, function(){$('.ad-hidden').show();});

void main()

{float x_y[N][2],A[N][K+1],B[N][M+1],AA[K+1],BB[M+1],a[M+1][M+2],m[M+1]; int i,j,n; printf("Please inqut %d\n",N); for(i=0;i

printf("(x%d y%d):",i,i);

scanf("%f %f",&x_y[i][0],&x_y[i][1]); }

for(i=0;i

A[i][0]=1;

for(j=1;j<=K;j++)

A[i][j]=A[i][j-1]*x_y[i][0]; for(j=0;j<=M;j++)

B[i][j]=A[i][j]*x_y[i][1]; }

for(j=0;j<=K;j++)

for(AA[j]=0,i=0;i

for(BB[j]=0,i=0;i

printf("系数矩阵是\n"); for(i=0;i<=n;i++) {for(j=0;j<=n+1;j++) printf("%f ",a[i][j]); printf("\n");} xiaoyuan(n,a); huidai(n,a,m);

printf("曲线方程是:\ny(x)=%g",m[0]); for(i=1;i<=n;i++) {

printf(" + %g",m[i]); for(j=0;j

printf("*X"); } }

}