数值分析多种算法C语言代码－推荐

1、离散傅立叶变换与反变换

/************************************************************************ 

* 离散傅立叶变换与反变换 

* 输入: 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,q,w,s; 

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]; 

} 

} 

}

 

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

/************************************************************************ 

* 用四阶亚当姆斯预估计求解初值问题，其中一阶微分方程未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); 

} 

3、几种常见随机数的产生

#include "stdlib.h" 

#include "stdio.h" 

#include "math.h" 

double uniform(double a,double b,long int* seed); 

double gauss(double mean,double sigma,long int *seed); 

double exponent(double beta,long int *seed); 

double laplace(double beta,long int* seed); 

double rayleigh(double sigma,long int *seed); 

double weibull(double a,double b,long int*seed); 

int bn(double p,long int*seed); 

int bin(int n,double p,long int*seed); 

int poisson(double lambda,long int *seed); 

void main() 

{ 

double a,b,x,mean; 

int i,j; 

long int s; 

a=4; 

b=0.7; 

s=13579; 

mean=0; 

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

{ 

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

{ 

x=poisson(a,&s); 

mean+=x; 

printf("%-13.7f",x); 

} 

printf("\n"); 

} 

mean/=50; 

printf("平均值为:%-13.7f\n",mean); 

} 

/******************************************************************* 

* 求[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); 

} 

/******************************************************************* 

* 正态分布 

* 输入: 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); 

} 

/******************************************************************* 

* 指数分布 

* 输入: 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); 

} 

/******************************************************************* 

* 拉普拉斯随机分布 

* 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); 

} 

/******************************************************************** 

* 瑞利分布 

* 

********************************************************************/ 

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); 

} 

/************************************************************************/ 

/* 韦伯分布                                                                     */ 

/************************************************************************/ 

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); 

} 

/************************************************************************/ 

/* 贝努利分布                                                           */ 

/************************************************************************/ 

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); 

} 

/************************************************************************/ 

/* 二项式分布                                                           */ 

/************************************************************************/ 

int bin(int n,double p,long int*seed) 

{ 

int i,x; 

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

x+=bn(p,seed); 

return(x); 

} 

/************************************************************************/ 

/* 泊松分布                                                             */ 

/************************************************************************/ 

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); 

} 

4、指数平滑法预测数据

/************************************************************************ 

* 本算法用指数平滑法预测数据 

* 输入: 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&
1978f
lt;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; 

} 

} 

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

精度比欧拉方法高 

但是感觉依然不理想 

/************************************************************************  * 用四阶(定步长)龙格--库塔法求解初值问题，其中一阶微分方程未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);  }

6、改进的欧拉方法求解微分方程初值问题

感觉精度比较低 

/************************************************************************  * 用改进的欧拉方法求解初值问题，其中一阶微分方程未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);  }

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

 

/************************************************************************  * 中心差分(矩形）公式计算函数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);  }

8。高斯10点法求积分

code] 

/******************************************************************** 

* 用高斯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))); 

} 

[/code] 

9、龙贝格法求积分

 

/********************************************************************  * 用龙贝格法计算函数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、复合辛普生法求积分

 

#include "stdio.h"  double composite_simpson(double(*f)(double),double a,double b,int n);  double myfun(double x);  void main()  {  double(*fun)(double);  double a,b,S4;  int n;  a=0;  b=1;  n=4;  fun=myfun;  S4=composite_simpson(fun,a,b,n);  printf("\n积分值为:%f\n",S4);  }  /********************************************************************  * 用复合辛普生法计算函数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);  }  double myfun(double x)  {  double y;  y=4.0/(1+x*x);  return(y);  }

11、最小二乘法拟合

 

/***************************************************************  * 本算法用最小二乘法依据指定的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 ,  double y ,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++)   /*选择第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);  }

12.埃特金插值法

/******************************************************  * 用埃特金插值法依据N个已知数据点计算函数值  * 输入: n--已知数据点的个数N-1  *       x--已知数据点第一坐标的N维列向量  * y--已知数据点第二坐标的N维列向量  * xx-插值点第一坐标  *       eps--求解精度  * 输出: 函数返回值所求插值点的第二坐标  ******************************************************/  double aitken(int n,double x ,double y ,double xx,double eps)  {  double d ;  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);  }  }

 

13、牛顿插值法

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

14、拉格朗日插值法

/******************************************************  * 用拉格朗日插值法依据N个已知数据点即使函数值  * 输入: n--已知数据点的个数N-1  *       x--已知数据点第一坐标的N维列向量  * y--已知数据点第二坐标的N维列向量  * xx-插值点第一坐标  * 输出: 函数返回值所求插值点的第二坐标  ******************************************************/  double lagrange(int n,double x ,double y ,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);  }

15、逆矩阵法求解矩阵方程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 ,int m,double b [M])  {  double det,f ;  int i,j,k;  det = gaussian_jodan(n,a);  if(det==0) return (0);  for(k=0;k<m;k++) /*做矩阵乘法AB*/  {  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 )  {  int i,j,k,mk;  int p ;  /*记录主行元素在原矩阵中的位置*/  double det,m,f;  det = 1.0;  for(k=0;k<n;k++)  {  m=a[k][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++) /*将第K列主元素换行到主对角线上*/  {  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);  }

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

/*************************************************************************  * 高斯列主元素消去法求解矩阵方程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 ,int m,double b [M])  {  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++)   /*选择第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);  }

17、任意多边形的面积计算（包括凹多边形的）

任意多边形的面积计算（包括凹多边形的，以及画多边形时线相交了的判别），最好提供相关资料，越详细越好，谢谢   

---------------------------------------------------------------   

我们都知道已知A（x1,y1）B（x2,y2）C（x3,y3）三点的面积公式为   

                      ¦x1    x2    x3  ¦   

S(A,B,C)  =       ¦y1    y2    y3  ¦  *  0.5    (当三点为逆时针时为正，顺时针则为负的)   

                      ¦1      1      1  ¦   

对多边形A1A2A3、、、An（顺或逆时针都可以），设平面上有任意的一点P，则有：   

  S（A1，A2，A3，、、、，An）   

  =  abs(S（P，A1，A2）  +  S（P，A2，A3）+、、、+S（P，An，A1）)   

P是可以取任意的一点，用（0，0）就可以了。   

这种方法对凸和凹多边形都适用。   

还有一个方法：   

任意一个简单多边形，当它的各个顶点位于网格的结点上时，它的面积数S=b/2+c+1   

其中：b代表该多边形边界上的网络结点数目   

          c代表该多边形内的网络结点数目  

所以把整个图形以象素为单位可以把整个图形分成若干个部分，计算该图形边界上的点b和内部的点c就得到面积数S了，然后把S乘以一个象素的面积就是所求的面积了。 

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

      非常感谢“山城棒棒儿的MATLAB&FPGA世界 ”，从中获益匪浅，有时间也会整理一些相关的数值分析的代码，共享给急切需要，没有目标的网友同行们！希望能在您能获多获少都会有所收获，那将是我最幸福的事！

转自 http://www.cnblogs.com/haohello/archive/2007/04/26/727744.html