相关数值分析多种算法代码

整理一些相关的数值分析的代码，共享给急切需要同行们！希望能在您能获多获少都会有所收获.>_<呵呵.

离散傅立叶变换与反变换

//************************************************************************// 离散傅立叶变换与反变换// 输入: x--要变换的数据的实部, y--要变换的数据的虚部//       a--变换结果的实部, b--变换结果的虚部//       n--数据长度//    sign--sign=1时，计算离散傅立叶正变换;sign=-1时;计算离散傅立叶反变换//************************************************************************void dft(double x[],double y[],double a[],double b[],int n,int sign){  int i,k;  double c,d,w,s;  double q= 6.28318530718/n;  for(k=0;k<n;k++){     w=k*q;     a[k]=b[k]=0.0;     for(i=0;i<n;i++){       d=i*w;       c=cos(d);       s=sin(d)*sign;       a[k]+=c*x+s*y;       b[k]+=c*y-s*x;     }  }  if(sign==-1){    c=1.0/n;    for(k=0;k<n;k++){      a[k]=c*a[k];      b[k]=c*b[k];    }  }} 

常用几种随机数产生的分布

泊松分布

int poisson(double lambda,long int *seed){  int i,x;  double a,b,u;    a=exp(-lambda);  i=0;  b=1.0;   do {    u=uniform(0.0,1.0,seed);    b*=u;    i++;  } while(b>=a);  x=i-1;  return(x);} 

贝努利分布

int bn(double p,long int*seed){  int x;  double u;  u=uniform(0.0,1.0,seed);  x=(u<=p)?1:0;  return(x);} 

韦伯分布

double weibull(double a,double b,long int*seed){  double u,x;    u=uniform(0.0,1.0,seed);  u=-log(u);  x=b*pow(u,1.0/a);    return(x);} 

瑞利分布

double rayleigh(double sigma,long int *seed){  double u,x;  u=uniform(0.,1.,seed);   x=-2.0*log(u);  x=sigma*sqrt(x);  return(x);} 

拉普拉斯随机分布

//******************************************************************// 拉普拉斯随机分布// beta--拉普拉斯分布的参数// *seed--随机数种子//******************************************************************double laplace(double beta,long int* seed){  double u1,u2,x;    u1 = uniform(0.,1.,seed);  u2 = uniform(0.,1.,seed);    if(u1<=0.5)     x = -beta*log(1.-u2);  else    x=beta*log(u2);    return(x);} 

指数分布

//******************************************************************// 指数分布// 输入: beta--指数分布均值//       seed--种子//******************************************************************double exponent(double beta,long int *seed){  double u,x;  u=uniform(0.0,1.0,seed);    x=-beta*log(u);  return(x);}

正态分布

//*******************************************************************// 正态分布// 输入: mean--双精度实型变量，正态分布的均值//      sigma--双精度实型变量，正态分布的均方差//       seed--长整型指针变量，*seed为随机数的种子  //*******************************************************************double gauss(double mean,double sigma,long int*seed){  int i;  double x,y;  for(x=0,i=0;i<12;i++)    x+=uniform(0.0,1.0,seed);  x=x-6.0;  y=mean+x*sigma;    return(y);}

均匀分布

//*******************************************************************// 求[a,b]上的均匀分布// 输入: a--双精度实型变量，给出区间的下限//       b--双精度实型变量，给出区间的上限//    seed--长整型指针变量，*seed为随机数的种子  //*******************************************************************double uniform(double a,double b,long int*seed){   double t; *seed=2045*(*seed)+1; *seed=*seed-(*seed/1048576)*1048576; t=(*seed)/1048576.0; t=a+(b-a)*t; return(t);}

几种求解初值问题的方法

四阶亚当姆斯预估计求解初值问题

//************************************************************************// 用四阶亚当姆斯预估计求解初值问题，其中一阶微分方程未y'=f(x,y)// 初始条件为x=x[0]时,y＝y[0].// 输入: f--函数f(x,y)的指针//       x--自变量离散值数组(其中x[0]为初始条件)//       y--对应于自变量离散值的函数值数组(其中y[0]为初始条件)//       h--计算步长//       n--步数// 输出: x为说求解的自变量离散值数组//       y为所求解对应于自变量离散值的函数值数组//************************************************************************double adams(double(*f)(double,double),double x[], double y[],double h,int n){  double dy[4],c,p,c1,p1,m;  int i,j;  runge_kuta(f,x,y,h,3);  for(i=0;i<4;i++)   dy=(*f)(x,y);  c=0.0;   p=0.0;  for(i=4;i<n+1;i++){    x=x[i-1]+h;    p1=y[i-1]+h*(55*dy[3]-59*dy[2]+37*dy[1]-9*dy[0])/24;    m=p1+251*(c-p)/270;    c1=y[i-1]+h*(9*(*f)(x,m)+19*dy[3]-5*dy[2]+dy[1])/24;    y=c1-19*(c1-p1)/270;    c=c1; p=p1;    for(j=0;j<3;j++)      dy[j]=dy[j+1];    dy[3]=(*f)(x,y);  }  return(0);} 

四阶(定步长)龙格--库塔法求解初值问题

/************************************************************************* 用四阶(定步长)龙格--库塔法求解初值问题，其中一阶微分方程未y'=f(x,y)* 初始条件为x=x[0]时,y＝y[0].* 输入: f--函数f(x,y)的指针*       x--自变量离散值数组(其中x[0]为初始条件)*       y--对应于自变量离散值的函数值数组(其中y[0]为初始条件)*       h--计算步长*       n--步数* 输出: x为说求解的自变量离散值数组*       y为所求解对应于自变量离散值的函数值数组************************************************************************/double runge_kuta(double(*f)(double,double),double x[], double y[],double h,int n){  int i;  double xs,ys,xp,yp,dy;  xs=x[0]+n*h;  for(i=0;i<n;i++){    ys=y;    dy=(*f)(x,y); //k1    y[i+1]=y+h*dy/6;    xp=x+h/2;    yp=ys+h*dy/2;    dy=(*f)(xp,yp); //k2    y[i+1]+=h*dy/3;    yp=ys+h*dy/2;    dy=(*f)(xp,yp);  //k3    y[i+1]+=h*dy/3;    xp+=h/2;    yp=ys+h*dy;    dy=(*f)(xp,yp); //k4    y[i+1]+=h*dy/6;    x[i+1]=xp;    if(x[i+1]>=xs)     return (0);  }  return(0);} 

用改进的欧拉方法求解初值问题

/************************************************************************* 用改进的欧拉方法求解初值问题，其中一阶微分方程未y'=f(x,y)* 初始条件为x=x[0]时,y＝y[0].* 输入: f--函数f(x,y)的指针*       x--自变量离散值数组(其中x[0]为初始条件)*       y--对应于自变量离散值的函数值数组(其中y[0]为初始条件)*       h--计算步长*       n--步数* 输出: x为说求解的自变量离散值数组*       y为所求解对应于自变量离散值的函数值数组************************************************************************/double proved_euler(double(*f)(double,double),double x[],double y[],double h,int n){  int i;  double xs,ys,yp;  for(i=0;i<n;i++){   ys=y;   xs=x;   y[i+1]=y;      yp=(*f)(xs,ys); //k1   y[i+1]+=yp*h/2.0;   ys+=h*yp;   xs+=h;      yp=(*f)(xs,ys); //k2   y[i+1]+=yp*h/2.0;   x[i+1]=xs;  }  return(0);}

中心差分(矩形）公式求导

/************************************************************************* 中心差分(矩形）公式计算函数f(x)在a点的导数值* 输入: f--函数f(x)的指针*       a--求导点*       h--初始步长*       eps--计算精度*       max_it--最大循环次数* 输出: 返回值为f(x)在a点的导数************************************************************************/double central_difference(double (*f)(double),double a, double h,double eps,int max_it){   double ff,gg;   int k;   ff=0.0;   for(k=0;k<max_it;k++){     gg=((*f)(a+h)-(*f)(a-h))/(h+h);     if(fabs(gg-ff)<eps)       return(gg);     h*=0.5;     ff=gg;   }   if(k==max_it){     printf("未能达到精度要求,需增大迭代次数!");     return(0);   }   return(gg);}

几种求积分的方法

龙贝格法求积分

/********************************************************************* 用龙贝格法计算函数f(x)从a到b的积分值* 输入: f--函数f(x)的指针*       a--积分下限*       b--积分上限*       eps--计算精度*       max_it--最大迭代次数* 输出: 返回值为f(x)从a点到b点的积分值*******************************************************************/double romberg(double(*f)(double),double a,double b,double eps,int max_it){  double *t,h;  int i,m,k;  if(!(t=(double *)malloc(max_it*sizeof(double)+1)))     return(ERROR_CODE);  h=b-a;  t[1]=h*((*f)(a)+(*f)(b))/2.0;  printf("%18.10e\n",t[1]);   for(k=2;k<max_it+1;k++){    double s,sm;    h*=0.5;    s=0.0;    for(i=0;i<pow(2,k-2);i++)       s+=(*f)(a+(2*i+1)*h);    sm=t[1];    t[1]=0.5*t[1]+h*s;    for(m=2;m<k+1;m++){      s=t[m];      t[m]=t[m-1]+(t[m-1]-sm)/(pow(4,m-1)-1);    if(m<k)     sm=s;  }  for(m=1;m<k+1;m++)    printf("%18.10e",t[m]);  printf("\n");  if(fabs(t[k]-sm)<eps){    sm=t[k];    free(t);    return(sm);  } } return(ERROR_CODE);}

高斯10点法求积分

/********************************************************************* 用高斯10点法计算函数f(x)从a到b的积分值* 输入: f--函数f(x)的指针*       a--积分下限*       b--积分上限* 输出: 返回值为f(x)从a点到b点的积分值*******************************************************************/double gauss_legendre(double(*f)(double),double a,double b){  const int n=10;  const double z[10]={-0.9739065285,-0.8650633677,-0.6794095683,                      -0.4333953941,-0.1488743390,0.1488743390,                      0.4333953941,0.6794095683,0.8650633677,                      0.9739065285};  const double w[10]={0.0666713443,0.1494513492,0.2190863625,                      0.2692667193,0.2955242247,0.2955242247,                      0.2692667193,0.2190863625,0.1494513492,                      0.0666713443};  double y,gg;  int i;  gg=0.0;  for(i=0;i<n;i++){    y=(z[i]*(b-a)+a+b)/2.0;    gg+=w[i]*(*f)((double)y);  }  return((double)((gg*(b-a)/2.0)));} 

复合辛普生法求积分

/********************************************************************* 用复合辛普生法计算函数f(x)从a到b的积分值* 输入: f--函数f(x)的指针*       a--积分下限*       b--积分上限*       n--分段数* 输出: 返回值为f(x)从a点到b点的积分值*******************************************************************/double composite_simpson(double(*f)(double),double a,double b,int n){  double s,h,x;  int i;  printf("x\t\tf(x)\t\ts\n");   s=(*f)(a)-(*f)(b);  h=(b-a)/(2*n);  x=a;  for(i=1;i<2*n;i+=2){    x+=h;    s+=4*(*f)(x);    printf("%f\t%f\t%f\n",x,(*f)(x),s*h/3);    x+=h;    s+=2*(*f)(x);    printf("%f\t%f\t%f\n",x,(*f)(x),s*h/3);  }  return(s*h/3);}

最小二乘法拟合

/**************************************************************** 本算法用最小二乘法依据指定的M个基函数及N个已知数据进行曲线拟和* 输入: m--已知数据点的个数M*       f--M维基函数向量* n--已知数据点的个数N-1*       x--已知数据点第一坐标的N维列向量*       y--已知数据点第二坐标的N维列向量*       a--无用* 输出: 函数返回值为曲线拟和的均方误差*       a为用基函数进行曲线拟和的系数,*       即a[0]f[0]+a[1]f[1]+...+a[M]f[M].****************************************************************/double mini_product(int m,double(*f[M])(double),int n,double x[N],                    double y[N],double a[M]){   double e,ff,b[M][M],c[M][1];   int i,j,k;   for(j=0;j<m;j++){    //计算最小均方逼近矩阵及常向量    for(k=0;k<m;k++){      b[j][k]=0.0;      for(i=0;i<n;i++)        b[j][k]+=(*f[j])(x)*(*f[k])(x);    }    c[j][0]=0.0;    for(i=0;i<n;i++)     c[j][0]+=(*f[j])(x)*y;   }   gaussian_elimination(m,b,1,c);   //求拟和系数   for(i=0;i<m;i++)     a=c[0];   e=0.0;   for(i=0;i<n;i++) { //计算均方误差     ff=0.0;     for(j=0;j<m;j++)       ff+=a[j]*(*f[j])(x);     e+=(y-ff)*(y-ff);   }  return(e);}

/************************************************************************** 高斯列主元素消去法求解矩阵方程AX=B,其中A是N*N的矩阵,B是N*M矩阵* 输入: n----方阵A的行数*       a----矩阵A*       m----矩阵B的列数*       b----矩阵B* 输出: det----矩阵A的行列式值*       a----A消元后的上三角矩阵*       b----矩阵方程的解X**************************************************************************/double gaussian_elimination(int n,double a[M][M],int m,double b[M][1]){  int i,j,k,mk;  double det,mm,f;  det = 1.0;   for(k = 0;k<n-1;k++){   mm=a[k][k];   mk = k;   for(i=k+1;i<n;i++) {    if(fabs(mm)<fabs(a[k])){     mm = a[k];     mk = i;    }   }  if(fabs(mm)<EPS)   return(0);  if(mk!=k){    for(j=k;j<n;j++){     f = a[k][j];     a[k][j]=a[mk][j];     a[mk][j]=f;    }    for(j=0;j<m;j++){     f = b[k][j];     b[k][j]=b[mk][j];     b[mk][j]=f;    }    det = -det;  }   for(i=k+1;i<n;i++){    mm = a[k]/a[k][k];    a[k]=0.0;    for(j=k+1;j<n;j++)      a[j]=a[j]-mm*a[k][j];    for(j=0;j<m;j++)      b[j]=b[j]-mm*b[k][j];  }  det = det*a[k][k]; } if(fabs(a[k][k])<EPS)   return 0; det=det*a[k][k]; for(i=0;i<m;i++) {   b[n-1]=b[n-1]/a[n-1][n-1];   for(j=n-2;j>=0;j--){     for(k=j+1;k<n;k++)       b[j]=b[j]-a[j][k]*b[k];   b[j]=b[j]/a[j][j];  } } return(det);}

几种插值法的方法

埃特金插值法

/******************************************************* 用埃特金插值法依据N个已知数据点计算函数值* 输入: n--已知数据点的个数N-1*       x--已知数据点第一坐标的N维列向量* y--已知数据点第二坐标的N维列向量* xx-插值点第一坐标*       eps--求解精度* 输出: 函数返回值所求插值点的第二坐标******************************************************/double aitken(int n,double x[N],double y[N],double xx,double eps){  double d[N];  int i,j;  for(i=0;i<=n;i++)    d=y;  for(i=0;i<=n;i++){    for(j=0;j<i;j++)     d=(d*(x[j]-xx)-d[j]*(x-xx))/(x[j]-x);    if(d-d[i-1]<eps)     return(d);  }}

牛顿插值法

/******************************************************* 用牛顿插值法依据N个已知数据点即使函数值* 输入: n--已知数据点的个数N-1*       x--已知数据点第一坐标的N维列向量* y--已知数据点第二坐标的N维列向量* xx-插值点第一坐标* 输出: 函数返回值所求插值点的第二坐标******************************************************/double newton(int n,double x[N],double y[N],double xx){  double d[N],b;  int i,j;  for(i=0;i<=n;i++)    d=y;  for(i=n-1;i>=0;i--)     for(j=i+1;j<=n;j++){      if(fabs(x-x[j])<EPS)       return 0;      d[j]=(d[j-1]-d[j])/(x-x[j]);    }  b=d[n];  for(i=n-1;i>=0;i--)    b=d+(xx-x)*b;  return b;} 

拉格朗日插值法

/******************************************************* 用拉格朗日插值法依据N个已知数据点即使函数值* 输入: n--已知数据点的个数N-1*       x--已知数据点第一坐标的N维列向量* y--已知数据点第二坐标的N维列向量* xx-插值点第一坐标* 输出: 函数返回值所求插值点的第二坐标******************************************************/double lagrange(int n,double x[N],double y[N],double xx){  double p,yy;  int i,j;  yy = 0.0;  for(i=0;i<=n;i++){   p=1.0;   for(j=0;j<=n;j++)    if(i!=j){     if(fabs(x-x[j])<EPS)       return 0;     p=p*(xx-x[j])/(x-x[j]);   }   yy=yy+p*y;  }  return(yy);}

求解矩阵方程AX=B的方法

逆矩阵法求解矩阵方程

/************************************************************************** 逆矩阵法求解矩阵方程AX=B，其中A是N*N的矩阵,B是N*N矩阵* 输入: n----方阵A的行数*       a----矩阵A*       m----矩阵B的列数*       b----矩阵B* 输出: det----矩阵A的行列式值*       a----A的逆矩阵*       b----矩阵方程的解X**************************************************************************/double gaussian_jodan_solve(int n,double a[N][N],int m,double b[N][M]){  double det,f[N];  int i,j,k;  det = gaussian_jodan(n,a);  if(det==0) return (0);  for(k=0;k<m;k++){   for(i=0;i<n;i++){    f=0.0;    for(j=0;j<n;j++)     f=f+a[j]*b[j][k];   }   for(i=0;i<n;i++)    b[k]=f;  }  return(det);} 

调用到的求逆矩阵的子函数

/************************************************************************** 高斯--约当列主元素法求矩阵方程A的逆矩阵，其中A是N*N的矩阵* 输入: n----方阵A的行数*       a----矩阵A* 输出: det--A的行列式的值*       a----A的逆矩阵**************************************************************************/double gaussian_jodan(int n,double a[N][N]){  int i,j,k,mk;  int p[N];   double det,m,f;  det = 1.0;  for(k=0;k<n;k++){    m=a[k][k];      mk=k;    for(i=k+1;i<n;i++)     if(fabs(m)<fabs(a[k])){        m=a[k];        mk=i;     }    if(fabs(m)<EPS) return(0);    if(mk!=k){       for(j=0;j<n;j++){        f=a[k][j];        a[k][j]=a[mk][j];        a[mk][j]=f;       }       p[k]=mk;       det = -det;    }    else      p[k]=k;    det=det*m;    for(j=0;j<n;j++)       if(j!=k)        a[k][j]=a[k][j]/a[k][k];    a[k][k]=1.0/a[k][k];    for(i=0;i<n;i++) {       if(i!=k){        for(j=0;j<n;j++)          if(j!=k)            a[j]=a[j]-a[k]*a[k][j];        a[k]=-a[k]*a[k][k];       }    }  }  for(k=n-2;k>=0;k--){   if(p[k]!=k)    for(i=0;i<n;i++){     f=a[k];     a[k]=a[p[k]];     a[p[k]]=f;    }  }  return(det);} 

高斯列主元素消去法求解矩阵方程

//************************************************************************// 高斯列主元素消去法求解矩阵方程AX=B,其中A是N*N的矩阵,B是N*M矩阵// 输入: n----方阵A的行数//       a----矩阵A//       m----矩阵B的列数//       b----矩阵B// 输出: det----矩阵A的行列式值//      a----A消元后的上三角矩阵//       b----矩阵方程的解X//*************************************************************************double gaussian_elimination(int n,double a[N][N],int m,double b[N][M]){  int i,j,k,mk;  double mm,f;  double det = 1.0;    for(k = 0;k<n-1;k++){  //选主元并消元   mm=a[k][k];   mk = k;   for(i=k+1;i<n;i++){ //选择第K列主元素     if(fabs(mm)<fabs(a[k])){       mm = a[k];       mk = i;     }   }     if(fabs(mm)<EPS)     return(0);   if(mk!=k){ //将第K列主元素换行到对角线上     for(j=k;j<n;j++){      f = a[k][j];      a[k][j]=a[mk][j];       a[mk][j]=f;     }         for(j=0;j<m;j++){      f = b[k][j];      b[k][j]=b[mk][j];      b[mk][j]=f;     }    det = -det;   }      for(i=k+1;i<n;i++){ //将第K列对角线以下消元为零     mm = a[k]/a[k][k];     a[k]=0.0;     for(j=k+1;j<n;j++)       a[j]=a[j]-mm*a[k][j];     for(j=0;j<m;j++)       b[j]=b[j]-mm*b[k][j];   }   det = det*a[k][k]; } if(fabs(a[k][k])<EPS)   return 0; det=det*a[k][k]; for(i=0;i<m;i++) {//回代求解   b[n-1]=b[n-1]/a[n-1][n-1];   for(j=n-2;j>=0;j--){     for(k=j+1;k<n;k++)        b[j]=b[j]-a[j][k]*b[k];b[j]=b[j]/a[j][j];   } } return(det);}

指数平滑法预测数据

//***********************************************************************// 本算法用指数平滑法预测数据// 输入: k--平滑周期//       n--原始数据个数//       m--预测步数//       alfa--加权系数//       x--指向原始数据数组指针// 输出: s1--返回值为指向一次平滑结果数组指针//       s2--返回值为指向二次指数平滑结果数组指针//       s3--返回值为指向三次指数平滑结果数组指针//      xx--返回值为指向预测结果数组指针//***********************************************************************void phyc(int k,int n,int m,double alfa,double x[N_MAX], double s1[N_MAX],double s2[N_MAX],double s3[N_MAX],double xx[N_MAX]){  double a,b,c,beta;  int i;  s1[k-1]=0;  for(i=0;i<k;k++)    s1[k-1]+=x;  s1[k-1]/=k;   for(i=k;i<=n;i++)    s1=alfa*x+(1-alfa)*s1[i-1];  s2[2*k-2]=0;   for(i=k-1;i<2*k-1;i++)    s2[2*k-2]+=s1;  s2[2*k-2]/=k;  for(i=2*k-1;i<=n;i++)    s2=alfa*s1+(1-alfa)*s2[i-1];  s3[3*k-3]=0;  for(i=2*k-2;i<3*k-2;i++)    s3[3*k-3]+=s2;  s3[3*k-3]/=k;  for(i=3*k-2;i<=n;i++)    s3=alfa*s2+(1-alfa)*s3[i-1];  beta=alfa/(2*(1-alfa)*(1-alfa));  for(i=3*k-3;i<=n;i++){    a=3*s1-3*s2+s3;    b=beta*((6-5*alfa)*s1-2*(5-4*alfa)*s2+(4-3*alfa)*s3);    c=beta*alfa*(s1-2*s2+s3);    xx=a+b*m+c*m*m;  }}

关于Image Engineering & Computer Vision的更多讨论与交流，敬请关注本博和新浪微博songzi_tea.