数值分析多种算法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<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[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)三点的面积公式为  
                      &brvbar;x1    x2    x3  &brvbar;  
S(A,B,C)  =       &brvbar;y1    y2    y3  &brvbar;  *  0.5    (当三点为逆时针时为正,顺时针则为负的)  
                      &brvbar;1      1      1  &brvbar;  

对多边形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乘以一个象素的面积就是所求的面积了。  

评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值