康建
数值逼近C程序
2018-1-19 17:44
阅读:2638
标签:教学, 计算方法, C语言, C语言

C语言中数组、指针、函数都是初学者不易掌握和应用的内容,数值分析中的函数逼近问题,也是较抽象的内容,这些内容对理工科学生来说都是应该必学必会的。

C语言编写数值分析中的算法,即有利于学生巩固计算机语言及数值分析知识,也有利于学生提高分析解决问题的能力,从而深入、形象地理解这两门课程内容,掌握算法的精髓。

本文重点介绍了用C语言编写曲线拟合的最小二乘法程序,其中用到了“函数指针数组”。这个程序能使学生更好地理解基函数、最佳平方逼近、法方程等概念,同时程序中也涉及到线性方程组求解、文件读写、数组参数传递、绘图等。


最后也给出对同一组数据进行分段线性插值、拉格朗日插值、三次样条插值、三角插值等的图像,这样有利于比较各种数值逼近方法,促进对知识的掌握及对算法综合运用。

函数的最佳平方逼近及曲线拟合的最小二乘法原理在此不多阐述,简单来说就是由一类基函数的线性组合来逼近某一函数,或者说在简单的基函数张成的函数空间中找一函数,在最佳平方意义上逼近所研究函数。其原理在几何上也有直观的意义,即等价于“垂线段最短”,如图:




这也是求解“法方程”这一词的意义吧,法方程可由“多元函数求极值必要条件”导出,解超定线性方程组Ax=b也是这原理,可归结为解法方程A’Ax=A’b

程序中定义了几个基函数,为了程序的通用性,定义了“函数指针数组”,

即:double (*pF[3])(double x)={fi0,fi1,fi2};

这是我写这篇博文要强调的,这里定义的pF是一个“函数指针数组”,fi0,fi1…是已定义的基函数。

pF[]*紧邻,[]优先级高,首先说明pF是个数组,加*说明是“指针数组”,括号后又有(double x)说明是指向函数的,这函数该有一个double型的参数(不写x也可),函数的返回类型也是double型。

以下是程序,有点长,“一肩担尽古今愁”,这是期末复习时演示的一个程序,可以复习练习很多内容。

(也没写太多的注释,如真需要可另解释)

#include "stdafx.h"

#include <windows.h>

#include <conio.h>

#include <stdlib.h>

#include <math.h>

#include <malloc.h>

#define PI 3.1415926535

typedef struct{

 double R;

 double I;

}COMPLEX;

int LUp(double *A,double *x,double *b,intn)

{       inti,r,k;

  double s;

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

      A[i*n+0]=A[i*n+0]/A[0*n+0];

for(r=1;r<n;r++)

       {      for(i=r;i<n;i++)

                     {      s=0;

                            for(k=0;k<r;k++)

                                   s+=A[r*n+k]*A[k*n+i];

                            A[r*n+i]=A[r*n+i]-s;

                     }

               for(i=r+1;i<n;i++)

                     {      s=0;

                            for(k=0;k<r;k++)

                                   s+=A[i*n+k]*A[k*n+r];

                            A[i*n+r]=(A[i*n+r]-s)/A[r*n+r];

                     }

        s=0;

        for(k=0;k<=r-1;k++)

               s+=A[r*n+k]*b[k];

        b[r]=b[r]-s;

       }

  x[n-1]=b[n-1]/A[(n-1)*n+n-1];

  for(i=n-2;i>=0;i--)

       {s=0;

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

               s+=A[i*n+k]*x[k];

        x[i]=(b[i]-s)/A[i*n+i];

       }

return0;

}

void zhuiganfa(int N,double a[],doubleb[],double c[],double f[],double X[])

{ double *bt=new double[N];

 double *y=new double[N];

 int i;

 bt[0]=c[0]/b[0];

 for(i=1;i<=N-2;i++)

        bt[i]=c[i]/(b[i]-a[i]*bt[i-1]);

 y[0]=f[0]/b[0];

 for(i=1;i<=N-1;i++)

        y[i]=(f[i]-a[i]*y[i-1])/(b[i]-a[i]*bt[i-1]);

  X[N-1]=y[N-1];

  for(i=N-2;i>=0;i--)

       X[i]=y[i]-bt[i]*X[i+1];

   delete [] bt,y;

}

double spline(int n,double x[],doubley[],double bj0,double bjn ,double xx )

{ int j;double tf1,tf2;double S;

double *a=new double[n];

double *b=new double[n];

double *c=new double[n];

double *d=new double[n];

double *h=new double[n];

double *M=new double[n];

  for(j=0;j<=n-2;j++)

             h[j]=x[j+1]-x[j];

      for(j=1;j<=n-2;j++)

       a[j]=(h[j-1])/(h[j-1] +  h[j]);

     a[n-1]=1;

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

       b[j]=2;

     for(j=1;j<=n-2;j++)

       c[j]=(h[j])/(h[j-1] +h[j]);  

       c[0]=1;

       for(j=1;j<=n-2;j++)

       {   tf1= (y[j+1]-y[j])/(h[j]);  

           tf2= (y[j]-y[j-1])/(h[j-1]);

               d[j]=6*( tf1 - tf2  )/( h[j-1] + h[j] );

       }

           d[0]=6*((y[1]-y[0])/(h[0]) -bj0 )/(h[0]);

              d[n-1]=6*(bjn-(y[n-1]-y[n-2])/(h[n-2]))/(h[n-2]);

     zhuiganfa(n,a,b, c,d,M);

    for(j=0;j<=n-2;j++)

             {if(xx>=x[j]&&xx<x[j+1])

                   {S=M[j]*(x[j+1]-xx)*(x[j+1]-xx)*(x[j+1]-xx)/(6* (h[j])  )

                     +M[j+1]*(xx-x[j])*(xx-x[j])*(xx-x[j]) / (6* (h[j]) )

                           +(  y[j]- M[j]*  (h[j])*(h[j]) /6   ) *  ((x[j+1]-xx)  /(h[j]) )

                           +(  y[j+1]- M[j+1]*  (h[j])*(h[j]) /6   ) * ( (xx-x[j])  /(h[j]) );

                    }

             }

     delete [] a,b,c,d,h,M;

  return S;

}

double lagrange_line(double x[],doubley[],int n,double xx)

{

  double yy=0;

  int k;

             for(k=0;k<n-1;k++)

                   {if(xx>=x[k]&&xx<=x[k+1])

                          {yy=y[k]*(x[k+1]-xx)/(x[k+1]-x[k])+y[k+1]*(xx-x[k])/(x[k+1]-x[k]);

                            break;

                           }

                    }

 return yy;

}

double lagrange(double x[],double y[],intn,double xx)

{  double yy=0,Lk;

  int j,k;

     for(k=0;k<n;k++){

             Lk=1;

             for(j=0;j<n;j++){

                 if(j!=k)Lk*=(xx-x[j])/(x[k]-x[j]);

             }

      yy+=y[k]*Lk;

     }

     return yy;

}

int readdata(double x[],double y[])

{

 FILE *fp;

 if((fp=fopen("1228101331.txt","r"))==NULL)

        {printf("cannot open file\n");

        exit(0) ;}

 int i=0;

 while(!feof(fp))

 {   fscanf(fp,"%lf,%lf\n",&x[i],&y[i]);

       i++;

  }

 fclose(fp);

 return i;

}

int drawfunction(double a,double b,doubleya,double yb,double (*f)(double ),char *win,int red)

{//a,bx的区间,ya,yb为估计的y的区间

     getch();

     HWND hWnd=FindWindow(NULL,win);

  HDC hDC=GetDC(hWnd);

  int i,n=800,iX,iY;

   double wx0=50,wx1=600,wy0=50,wy1=500;

  double x,y;

     for(i=0;i<n;i++){

             x=a+i*(b-a)/n;

             y=f(a+i*(b-a)/n);

      iX=(int)(wx0+ (x-a)* (wx1-wx0)/(b-a) );

            iY=600-(int)(wy0+(y-ya)*(wy1-wy0)/(yb-ya) );

             SetPixel(hDC,iX,iY,RGB(red,0,255) );

     }

     getch();

     return 0;

}

double mySpline(double t)

{

     int N=16,n;

     double *x=new double [N];

     double *y=new double [N];

  double bj0=3.0,bjn=-0.1;

  n=readdata(x,y);

     if(t<x[0]||t>x[n-1])

             return 0;

     else

  {      bj0=(y[1]-y[0])/(x[1]-x[0]);

       bjn=(y[n-1]-y[n-2])/(x[n-1]-x[n-2]);

     return spline(n,x, y, bj0, bjn ,t);

  }

}

double myLagrange_line(double t)

{

     int N=16,n;

     double *x=new double [N];

     double *y=new double [N];

  n=readdata(x,y);

     if(t<x[0]||t>x[n-1])

             return 0;

     else

     return lagrange_line(x, y,n,t);

}

double mylagrange(double t)

{

     int N=16,n;

     double *x=new double [N];

     double *y=new double [N];

     n=readdata(x,y);

     if(t<x[0]||t>x[n-1])

             return 0;

     else

     return lagrange(x, y,n,t);

}

double fi0(double t){   return 1;}

double fi1(double t){   return t;}

double fi2(double t){   return t*t;}

double myLS_1p(double t)

{  int N=16,n;

  double *x=new double [N];

  double *y=new double [N];

  n=readdata(x,y);

  int p=2;

  double (*pF[2])(double t);

pF[0]=fi0;//指向基函数

pF[1]=fi1;

doubleG[2][2]={0},b[2]={0};

doublea[2];

inti,r,c;

for(i=0;i<n;i++)

{

       for(r=0;r<p;r++)

       {      for(c=0;c<p;c++)

                   G[r][c]+=pF[r](x[i])*pF[c](x[i]);

           b[r]+=pF[r](x[i])*y[i];

       }

}

 LUp(G[0],a,b,p);

doubles=0;

for(i=0;i<p;i++)

{   s+=a[i]*pF[i](t);

}

returns;

}

double myLS_2p(double t)

{

  int N=16,n;

  double *x=new double [N];

  double *y=new double [N];

  n=readdata(x,y);

  int p=3;

  double (*pF[3])(double x)={fi0,fi1,fi2};

doubleG[3][3]={0},b[3]={0};

doublea[3];

inti,r,c;

for(i=0;i<n;i++)

{

       for(r=0;r<p;r++)

       {      for(c=0;c<p;c++)

                   G[r][c]+=pF[r](x[i])*pF[c](x[i]);

           b[r]+=pF[r](x[i])*y[i];

       }

}

LUp(G[0],a,b,p);

doubles=0;

for(i=0;i<p;i++)

{   s+=a[i]*pF[i](t);

}

returns;

}

double f0(double t){   return 1;}

double f1(double t){   return cos(t);}

double f2(double t){   return sin(t);}

double f3(double t){   return cos(2*t);}

double f4(double t){   return sin(2*t);}

double f5(double t){   return cos(3*t);}

double f6(double t){   return sin(3*t);}

double f7(double t){   return cos(4*t);}

double myLS_sincos(double t)

{   int N=16,n;

  double *x=new double [N];

  double *y=new double [N];

  n=readdata(x,y);

  int p=8;

  double (*pF[8])(double x);

pF[0]=f0;

pF[1]=f1;    pF[2]=f2;

   pF[3]=f3;     pF[4]=f4;

   pF[5]=f5;     pF[6]=f6;

   pF[7]=f7;

doubleG[8][8]={0},b[8]={0};

doublea[8];

inti,r,c;

for(i=0;i<n;i++)

{

       for(r=0;r<p;r++)

       {      for(c=0;c<p;c++)

                   G[r][c]+=pF[r](x[i])*pF[c](x[i]);

           b[r]+=pF[r](x[i])*y[i];

       }

}

LUp(G[0],a,b,p);

doubles=0;

for(i=0;i<p;i++)

{   s+=a[i]*pF[i](t);

}

returns;

}

////fft

COMPLEX add(COMPLEX a,COMPLEX b)

{   COMPLEX c;

   c.R=a.R+b.R;     c.I=a.I+b.I;        return c;

}

COMPLEX sub(COMPLEX a,COMPLEX b)

{   COMPLEX c;

   c.R=a.R-b.R;     c.I=a.I-b.I; return c;

}

COMPLEX mul(COMPLEX a,COMPLEX b)

{   COMPLEX c;

   c.R=a.R*b.R-a.I*b.I;    c.I=a.R*b.I+b.R*a.I;

return c;

}

int nx(int k,int p)

{

  int q,t=0;

for(q=0;q<p;q++)  

              if((1<<q)&k)    

                     t+=1<<(p-q-1);

returnt;

}

void kj_FFT(COMPLEX *T,COMPLEX *F,int p)

{      

 intN,q,k,r,j,i,t;

COMPLEX *A1,*A2, *W;

double a;

  N=1<<p;

  A1=(COMPLEX *)malloc(sizeof(COMPLEX)*N);

  A2=(COMPLEX *)malloc(sizeof(COMPLEX)*N);

   W=(COMPLEX *)malloc(sizeof(COMPLEX)*N/2);

 for(i=0;i<N;i++)

    A1[i]=T[i];

 for(i=0;i<N/2;i++)

       {  a=-i*2*PI/N;  

          W[i].R=cos(a);  

          W[i].I=sin(a);

       }

 for(q=0;q<p;q++)

  {for(k=0;k<1<<q;k++)

          { r=1<<(p-q-1);

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

                     {t=k*( 1<<(p-q)  );

                       A2[j+t] =add(A1[j+t], mul(A1[j+t+r],W[ nx((j+t)>>(p-q-1),p)]) );

                       A2[j+t+r]=sub(A1[j+t], mul(A1[j+t+r],W[nx((j+t)>>(p-q-1),p)]) );

                      }

          }

       for(i=0;i<N;i++)

              A1[i]=A2[i];

  }

 for(i=0;i<N;i++)

 {  

  F[i].R=A2[nx(i,p)].R/N;

     F[i].I=A2[nx(i,p)].I/N;

  }

}

double myFFTsincos(double t)

{

     int N=16,n,p=4,j,k;

     double *x=new double [N];

     double *y=new double [N];

     n=readdata(x,y);

  COMPLEX *pT,*pF;

     pT=(COMPLEX *)malloc(sizeof(COMPLEX)*n);

     pF=(COMPLEX *)malloc(sizeof(COMPLEX)*n);

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

  {  pT[j].R=y[j];

          pT[j].I=0.;

  }

   kj_FFT(pT,pF,p);

              for(k=0;k<=n/2;k++)

              {   pF[k].R=pF[k].R*2;

                     pF[k].I=pF[k].I*2;

              }

     if(t<x[0]||t>x[n-1])

             return 0;

     else

     {      double s,tt;

              tt=2*PI*(t-x[0])/(x[n-1]-x[0]);

              s=pF[0].R/2;

              for(j=1;j<=n/2-1;j++)

              {s+=pF[j].R*cos(j*tt)-pF[j].I*sin(j*tt);

              }

                s=s+1/2*pF[n/2].R*cos(n/2*tt);

             return s;

  }

}

int main(int argc, char* argv[])

{     double xa=0,xb=3.75,ya=-0.2,yb=2;

 drawfunction(xa,xb,ya,yb,mySpline,argv[0],0);

     drawfunction(xa,xb,ya,yb,myLagrange_line,argv[0],50);

 drawfunction(xa,xb,ya,yb,myLS_1p,argv[0],150);

 drawfunction(xa,xb,ya,yb,myLS_2p,argv[0],200);

 drawfunction(xa,xb,ya,yb,myLS_sincos,argv[0],250);

 drawfunction(xa,xb,ya,yb,myFFTsincos,argv[0],0);

     return 0;

}

下图是两组数据的运行结果。






转载本文请联系原作者获取授权,同时请注明本文来自康建科学网博客。

链接地址:http://wap.sciencenet.cn/blog-797552-1095762.html?mobile=1

收藏

分享到:

当前推荐数:2
推荐人:
推荐到博客首页
网友评论2 条评论
确定删除指定的回复吗?
确定删除本博文吗?