曲线拟合

拉格朗日插值法:
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
 
doubleLAG(int,double *,double *,double);
 
voidmain()
{
         intn;
         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);
}
 
doubleLAG(n,x,y,t)
intn;
double *x;
double *y;
doublet;
{
         inti,j;
         doublep,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 <stdlib.h>
#include <malloc.h>
#include <math.h>
 
Smooth(double *,double *,double *,int,int,
            double *,double *,double *);
voidmain()
{
         inti,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。返回m1次拟合多项式的m个系数*/
                   intn;                                               /*整型变量,输入参数,给定数据点的个数*/
                   intm;                                               /*整型变量,输入参数,拟合多项式的项数*/
                   double *dt1;                            /*实型变量,输出参数,拟合多项式与数据点偏差的平方和*/
                   double *dt2;                            /*实型变量,输出参数,拟合多项式与数据点偏差的绝对值之和*/
                   double *dt3;                            /*实型变量,输出参数,拟合多项式与数据点偏差的绝对值最大值*/
         {
                   inti,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);
         }
 
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值