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

最新推荐文章于 2021-05-19 15:48:40 发布

Augusdi

最新推荐文章于 2021-05-19 15:48:40 发布

阅读量7.9k

点赞数 7

分类专栏： DSP数字信号处理

DSP数字信号处理专栏收录该内容

105 篇文章 23 订阅

订阅专栏

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<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点法求积分

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

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[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++)   /*选择第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[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); 
	} 
}

13、牛顿插值法

/****************************************************** 
* 用牛顿插值法依据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; 
}

14、拉格朗日插值法

/****************************************************** 
* 用拉格朗日插值法依据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); 
}

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[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++) /*做矩阵乘法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[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];  /*选第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[N][N],int m,double b[N][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乘以一个象素的面积就是所求的面积了。

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