是不是想寫成為y=at+bt^2+ct^3
但是實際擬合出來的表達式為y=a[3]+a[2]t+a[1]t^2+a[0]t^3會有個常數項的。
簡單的講,所謂擬合是指已知某函數的若幹離散函數值{f1,f2,…,fn},通過調整該函數中若幹待定系數f(λ1, λ2,…,λ3), 使得該函數與已知點集的差別(最小二乘意義)最小。如果待定函數是線性,就叫線性擬合或者線性回歸(主要在統計中),否則叫作非線性擬合或者非線性回歸。表達式也可以是分段函數,這種情況下叫作樣條擬合。
曲線擬合:
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
Smooth(double *x,double *y,double *a,int n,int m,double *dt1,double *dt2,double *dt3);
void main()
{
int i ,n ,m ;
double *x,*y,*a,dt1,dt2,dt3,b;
n = 12;// 12個樣點
m = 4; //3次多項式擬合
b = 0; //x的初值為0
/*分別為x,y,a分配存貯空間*/
x = (double *)calloc(n,sizeof(double));
if(x == NULL)
{
printf("內存分配失敗\n");
exit (0);
}
y = (double *)calloc(n,sizeof(double));
if(y == NULL)
{
printf("內存分配失敗\n");
exit (0);
}
a = (double *)calloc(n,sizeof(double));
if(a == NULL)
{
printf("內存分配失敗\n");
exit (0);
}
for(i=1;i<=n;i++)
{
x[i-1]=b+(i-1)*5;
/*每隔5取壹個點,這樣連續取12個點*/
}
y[0]=0;
y[1]=1.27;
y[2]=2.16;
y[3]=2.86;
y[4]=3.44;
y[5]=3.87;
y[6]=4.15;
y[7]=4.37;
y[8]=4.51;
y[9]=4.58;
y[10]=4.02;
y[11]=4.64;
/*x[i-1]點對應的y值是擬合已知值*/
Smooth(x,y,a,n,m,&dt1,&dt2,&dt3); /*調用擬合函數*/
for(i=1;i<=m;i++)
printf("a[%d] = %.10f\n",(i-1),a[i-1]);
printf("擬合多項式與數據點偏差的平方和為:\n");
printf("%.10e\n",dt1);
printf("擬合多項式與數據點偏差的絕對值之和為:\n");
printf("%.10e\n",dt2);
printf("擬合多項式與數據點偏差的絕對值最大值為:\n");
printf("%.10e\n",dt3);
free(x); /*釋放存儲空間*/
free(y); /*釋放存儲空間*/
free(a); /*釋放存儲空間*/
}
Smooth(double *x,double *y,double *a,int n,int m,double *dt1,double *dt2,double *dt3)//(x,y,a,n,m,dt1,dt2,dt3 )
//double *x; /*實型壹維數組,輸入參數,存放節點的xi值*/
//double *y; /*實型壹維數組,輸入參數,存放節點的yi值*/
//double *a; /*雙精度實型壹維數組,長度為m。返回m壹1次擬合多項式的m個系數*/
//int n; /*整型變量,輸入參數,給定數據點的個數*/
//int m; /*整型變量,輸入參數,擬合多項式的項數*/
//double *dt1; /*實型變量,輸出參數,擬合多項式與數據點偏差的平方和*/
//double *dt2; /*實型變量,輸出參數,擬合多項式與數據點偏差的絕對值之和*/
//double *dt3; /*實型變量,輸出參數,擬合多項式與數據點偏差的絕對值最大值*/
{
int i ,j ,k ;
double *s,*t,*b,z,d1,p,c,d2,g,q,dt;
/*分別為s ,t ,b分配存貯空間*/
s = (double *)calloc(n,sizeof(double));
if(s == NULL)
{
printf("內存分配失敗\n");
exit (0);
}
t = (double *)calloc(n,sizeof(double));
if(t == NULL)
{
printf("內存分配失敗\n");
exit (0);
}
b = (double *)calloc(n,sizeof(double));
if(b == NULL)
{
printf("內存分配失敗\n");
exit (0);
}
z = 0;
for(i=1;i<=n;i++)
z=z+x[i-1]/n; /*z為各個x的平均值*/
b[0]=1;
d1=n;
p=0;
c=0;
for(i=1;i<=n;i++)
{
p=p+x[i-1]-z;
c=c+y[i-1];
}
c=c/d1;
p=p/d1;
a[0]=c*b[0];
if(m>1)
{
t[1]=1;
t[0]=-p;
d2=0;
c=0;
g=0;
for(i=1;i<=n;i++)
{
q=x[i-1]-z-p;
d2=d2+q*q;
c=y[i-1]*q+c;
g=(x[i-1]-z)*q*q+g;
}
c=c/d2;
p=g/d2;
q=d2/d1;
d1=d2;
a[1]=c*t[1];
a[0]=c*t[0]+a[0];
}
for(j=3;j<=m;j++)
{
s[j-1]=t[j-2];
s[j-2]=-p*t[j-2]+t[j-3];
if(j>=4)
for(k=j-2;k>=2;k--)
s[k-1]=-p*t[k-1]+t[k-2]-q*b[k-1];
s[0]=-p*t[0]-q*b[0];
d2=0;
c=0;
g=0;
for(i=1;i<=n;i++)
{
q=s[j-1];
for(k=j-1;k>=1;k--)
q=q*(x[i-1]-z)+s[k-1];
d2=d2+q*q;
c=y[i-1]*q+c;
g=(x[i-1]-z)*q*q+g;
}
c=c/d2;
p=g/d2;
q=d2/d1;
d1=d2;
a[j-1]=c*s[j-1];
t[j-1]=s[j-1];
for(k=j-1;k>=1;k--)
{
a[k-1]=c*s[k-1]+a[k-1];
b[k-1]=t[k-1];
t[k-1]=s[k-1];
}
}
*dt1=0;
*dt2=0;
*dt3=0;
for(i=1;i<=n;i++)
{
q=a[m-1];
for(k=m-1;k>=1;k--)
q=q*(x[i-1]-z)+a[k-1];
dt=q-y[i-1];
if(fabs(dt)>*dt3)
*dt3=fabs(dt);
*dt1=*dt1+dt*dt;
*dt2=*dt2+fabs(dt);
}
/*釋放存儲空間*/
free(s);
free(t);
free(b);
return(1);
}