c.txt

用c语言写的数值分析的多种算法的代码。实现了多种数值分析算法。

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