兩個簡單的曲線擬合的方法：
拉格朗日插值法：
#include <stdio.h>
#include
#include <malloc.h>
double LAG(int,double *,double *,double);
void main()
{
int n;
double *x,*y,t,lag;
t = 0.15;
n = 6;
x = (double*)calloc(n,sizeof(double));
if(x == NULL)
{
printf("內存分配失敗/n");
exit(1);
}
y = (double*)calloc(n,sizeof(double));
if(y == NULL)
{
printf("內存分配失敗/n");
exit(1);
}
x[0] = 0;
x[1] = 0.1;
x[2] = 0.195;
x[3] = 0.3;
x[4] = 0.401;
x[5] = 0.5;
y[0] = 0.39894;
y[1] = 0.39695;
y[2] = 0.39142;
y[3] = 0.38138;
y[4] = 0.36812;
y[5] = 0.35206;
lag = LAG(n,x,y,t);
printf("拉各朗日插值後得到的結果是：/n");
printf("f(%.2f)=%e/n",t,lag);
free(x);
free(y);
}
double LAG(n,x,y,t)
int n;
double *x;
double *y;
double t;
{
int i,j;
double p,s;
s = 0;
for(i=0;i<n-1;i++)
{
p = 1;
for(j=0;j<n-1;j++)
if(i!=j)
p*=(t-x[j])/(x[i]-x[j]);
s+=p*y[i];
}
return (s);
}
曲線擬合：
#include <stdio.h>
#include
#include <malloc.h>
#include
Smooth(double *,double *,double *,int,int,
double *,double *,double *);
void main()
{
int i,n,m;
double *x,*y,*a,dt1,dt2,dt3,b;
n = 20;
m = 6;
b = 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)*0.1;
/*每隔0.1取一個點，這樣連續取n個點*/
y[i-1]=x[i-1]-exp(-x[i-1]);
/*計算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(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);
}
沒有留言:
張貼留言