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

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

分类: DSP数字信号处理 157人阅读 评论(0) 收藏 举报

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

  1. /************************************************************************  
  2. * 离散傅立叶变换与反变换  
  3. * 输入: x--要变换的数据的实部  
  4. * y--要变换的数据的虚部  
  5. *       a--变换结果的实部  
  6. *       b--变换结果的虚部  
  7. *       n--数据长度  
  8. *    sign--sign=1时,计算离散傅立叶正变换;sign=-1时;计算离散傅立叶反变换  
  9. ************************************************************************/   
  10. void dft(double x[], double y[], double a[], double b[], int n, int sign)   
  11. {   
  12.     int i,k;   
  13.     double c, d, q, w, s;   
  14.   
  15.     q=6.28318530718 / n;   
  16.     for(k = 0; k < n; k++)   
  17.     {   
  18.         w = k * q;   
  19.         a[k] = b[k] = 0.0;   
  20.         for(i = 0; i < n; i++)   
  21.         {   
  22.             d = i * w;   
  23.             c = cos(d);   
  24.             s = sin(d) * sign;   
  25.             a[k] += c * x + s * y;   
  26.             b[k] += c * y - s * x;   
  27.         }   
  28.     }   
  29.     if(sign == -1)   
  30.     {   
  31.         c = 1.0 / n;   
  32.         for(k = 0; k < n; k++)   
  33.         {   
  34.             a[k] = c * a[k];   
  35.             b[k] = c * b[k];   
  36.         }   
  37.     }   
  38. }   
/************************************************************************ 
* 离散傅立叶变换与反变换 
* 输入: 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四阶亚当姆斯预估计求解初值问题

  1. /************************************************************************  
  2. * 用四阶亚当姆斯预估计求解初值问题,其中一阶微分方程未y'=f(x,y)  
  3. * 初始条件为x=x[0]时,y=y[0].  
  4. * 输入: f--函数f(x,y)的指针  
  5. *       x--自变量离散值数组(其中x[0]为初始条件)  
  6. *       y--对应于自变量离散值的函数值数组(其中y[0]为初始条件)  
  7. *       h--计算步长  
  8. *       n--步数  
  9. * 输出: x为说求解的自变量离散值数组  
  10. *       y为所求解对应于自变量离散值的函数值数组  
  11. ************************************************************************/   
  12. double adams(double(*f)(double,double),double x[], double y[],double h,int n)   
  13. {   
  14.     double dy[4],c,p,c1,p1,m;   
  15.     int i,j;   
  16.   
  17.     runge_kuta(f,x,y,h,3);   
  18.     for(i=0;i<4;i++)   
  19.         dy=(*f)(x,y);   
  20.     c=0.0; p=0.0;   
  21.     for(i=4;i<n+1;i++)   
  22.     {   
  23.         x=x[i-1]+h;   
  24.         p1=y[i-1]+h*(55*dy[3]-59*dy[2]+37*dy[1]-9*dy[0])/24;   
  25.         m=p1+251*(c-p)/270;   
  26.         c1=y[i-1]+h*(9*(*f)(x,m)+19*dy[3]-5*dy[2]+dy[1])/24;   
  27.         y=c1-19*(c1-p1)/270;   
  28.         c=c1; p=p1;   
  29.         for(j=0;j<3;j++)   
  30.             dy[j]=dy[j+1];   
  31.         dy[3]=(*f)(x,y);   
  32.     }   
  33.     return(0);   
  34. }   
/************************************************************************ 
* 用四阶亚当姆斯预估计求解初值问题,其中一阶微分方程未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、几种常见随机数的产生

  1. #include "stdlib.h"    
  2. #include "stdio.h"    
  3. #include "math.h"    
  4.   
  5. double uniform(double a,double b,long int* seed);   
  6. double gauss(double mean,double sigma,long int *seed);   
  7. double exponent(double beta,long int *seed);   
  8. double laplace(double beta,long int* seed);   
  9. double rayleigh(double sigma,long int *seed);   
  10. double weibull(double a,double b,long int*seed);   
  11. int bn(double p,long int*seed);   
  12. int bin(int n,double p,long int*seed);   
  13. int poisson(double lambda,long int *seed);   
  14.   
  15. void main()   
  16. {   
  17.     double a,b,x,mean;   
  18.     int i,j;   
  19.     long int s;   
  20.   
  21.     a=4;   
  22.     b=0.7;   
  23.     s=13579;   
  24.     mean=0;   
  25.     for(i=0;i<10;i++)   
  26.     {   
  27.         for(j=0;j<5;j++)   
  28.         {   
  29.             x=poisson(a,&s);   
  30.             mean+=x;   
  31.             printf("%-13.7f",x);   
  32.         }   
  33.         printf("\n");   
  34.     }   
  35.     mean/=50;   
  36.     printf("平均值为:%-13.7f\n",mean);   
  37. }   
  38.   
  39. /*******************************************************************  
  40. * 求[a,b]上的均匀分布  
  41. * 输入: a--双精度实型变量,给出区间的下限  
  42. *       b--双精度实型变量,给出区间的上限  
  43. *    seed--长整型指针变量,*seed为随机数的种子    
  44. ********************************************************************/   
  45. double uniform(double a,double b,long int*seed)   
  46. {   
  47.     double t;   
  48.     *seed=2045*(*seed)+1;   
  49.     *seed=*seed-(*seed/1048576)*1048576;   
  50.     t=(*seed)/1048576.0;   
  51.     t=a+(b-a)*t;   
  52.   
  53.     return(t);   
  54. }   
  55.   
  56. /*******************************************************************  
  57. * 正态分布  
  58. * 输入: mean--双精度实型变量,正态分布的均值  
  59. *      sigma--双精度实型变量,正态分布的均方差  
  60. *       seed--长整型指针变量,*seed为随机数的种子    
  61. ********************************************************************/   
  62. double gauss(double mean,double sigma,long int*seed)   
  63. {   
  64.     int i;   
  65.     double x,y;   
  66.     for(x=0,i=0;i<12;i++)   
  67.         x+=uniform(0.0,1.0,seed);   
  68.     x=x-6.0;   
  69.     y=mean+x*sigma;   
  70.   
  71.     return(y);   
  72. }   
  73.   
  74. /*******************************************************************  
  75. * 指数分布  
  76. * 输入: beta--指数分布均值  
  77. *       seed--种子  
  78. *******************************************************************/   
  79. double exponent(double beta,long int *seed)   
  80. {   
  81.     double u,x;   
  82.     u=uniform(0.0,1.0,seed);   
  83.     x=-beta*log(u);   
  84.   
  85.     return(x);   
  86. }   
  87.   
  88. /*******************************************************************  
  89. * 拉普拉斯随机分布  
  90. * beta--拉普拉斯分布的参数  
  91. * *seed--随机数种子  
  92. *******************************************************************/   
  93. double laplace(double beta,long int* seed)   
  94. {   
  95.     double u1,u2,x;   
  96.   
  97.     u1=uniform(0.,1.,seed);   
  98.     u2=uniform(0.,1.,seed);   
  99.     if(u1<=0.5)   
  100.         x=-beta*log(1.-u2);   
  101.     else   
  102.         x=beta*log(u2);   
  103.   
  104.     return(x);   
  105. }   
  106.   
  107.   
  108.   
  109. /********************************************************************  
  110. * 瑞利分布  
  111.  
  112. ********************************************************************/   
  113. double rayleigh(double sigma,long int *seed)   
  114. {   
  115.     double u,x;   
  116.     u=uniform(0.,1.,seed);   
  117.     x=-2.0*log(u);   
  118.     x=sigma*sqrt(x);   
  119.     return(x);   
  120. }   
  121.   
  122. /************************************************************************/   
  123. /* 韦伯分布                                                                     */   
  124. /************************************************************************/   
  125. double weibull(double a,double b,long int*seed)   
  126. {   
  127.     double u,x;   
  128.   
  129.     u=uniform(0.0,1.0,seed);   
  130.     u=-log(u);   
  131.     x=b*pow(u,1.0/a);   
  132.   
  133.     return(x);   
  134. }   
  135.   
  136. /************************************************************************/   
  137. /* 贝努利分布                                                           */   
  138. /************************************************************************/   
  139. int bn(double p,long int*seed)   
  140. {   
  141.     int x;   
  142.     double u;   
  143.     u=uniform(0.0,1.0,seed);   
  144.     x=(u<=p)?1:0;   
  145.     return(x);   
  146. }   
  147.   
  148. /************************************************************************/   
  149. /* 二项式分布                                                           */   
  150. /************************************************************************/   
  151. int bin(int n,double p,long int*seed)   
  152. {   
  153.     int i,x;   
  154.     for(x=0,i=0;i<n;i++)   
  155.     x+=bn(p,seed);   
  156.     return(x);   
  157. }   
  158.   
  159. /************************************************************************/   
  160. /* 泊松分布                                                             */   
  161. /************************************************************************/   
  162. int poisson(double lambda,long int *seed)   
  163. {   
  164.     int i,x;   
  165.     double a,b,u;   
  166.   
  167.     a=exp(-lambda);   
  168.     i=0;   
  169.     b=1.0;   
  170.     do {   
  171.         u=uniform(0.0,1.0,seed);   
  172.         b*=u;   
  173.         i++;   
  174.     } while(b>=a);   
  175.     x=i-1;   
  176.     return(x);   
  177. }   
#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、指数平滑法预测数据

  1. /************************************************************************  
  2. * 本算法用指数平滑法预测数据  
  3. * 输入: k--平滑周期  
  4. *       n--原始数据个数  
  5. *       m--预测步数  
  6. *       alfa--加权系数  
  7. *       x--指向原始数据数组指针  
  8. * 输出: s1--返回值为指向一次平滑结果数组指针  
  9. *       s2--返回值为指向二次指数平滑结果数组指针  
  10. *       s3--返回值为指向三次指数平滑结果数组指针  
  11. *       xx--返回值为指向预测结果数组指针  
  12. ************************************************************************/   
  13. void phyc(int k,int n,int m,double alfa,double x[N_MAX],   
  14.  double s1[N_MAX],double s2[N_MAX],double s3[N_MAX],double xx[N_MAX])   
  15. {   
  16.     double a,b,c,beta;   
  17.     int i;   
  18.   
  19.     s1[k-1]=0;   
  20.     for(i=0;i<k;k++)   
  21.         s1[k-1]+=x;   
  22.     s1[k-1]/=k;   
  23.     for(i=k;i<=n;i++)   
  24.         s1=alfa*x+(1-alfa)*s1[i-1];   
  25.     s2[2*k-2]=0;   
  26.     for(i=k-1;i<2*k-1;i++)   
  27.         s2[2*k-2]+=s1;   
  28.     s2[2*k-2]/=k;   
  29.     for(i=2*k-1;i<=n;i++)   
  30.         s2=alfa*s1+(1-alfa)*s2[i-1];   
  31.     s3[3*k-3]=0;   
  32.     for(i=2*k-2;i<3*k-2;i++)   
  33.         s3[3*k-3]+=s2;   
  34.     s3[3*k-3]/=k;   
  35.     for(i=3*k-2;i<=n;i++)   
  36.         s3=alfa*s2+(1-alfa)*s3[i-1];   
  37.     beta=alfa/(2*(1-alfa)*(1-alfa));   
  38.     for(i=3*k-3;i<=n;i++)   
  39.     {   
  40.         a=3*s1-3*s2+s3;   
  41.         b=beta*((6-5*alfa)*s1-2*(5-4*alfa)*s2+(4-3*alfa)*s3);   
  42.         c=beta*alfa*(s1-2*s2+s3);   
  43.         xx=a+b*m+c*m*m;   
  44.     }   
  45. }   
/************************************************************************ 
* 本算法用指数平滑法预测数据 
* 输入: 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、四阶(定步长)龙格--库塔法求解微分初值问题

精度比欧拉方法高
但是感觉依然不理想

  1. /************************************************************************  
  2. * 用四阶(定步长)龙格--库塔法求解初值问题,其中一阶微分方程未y'=f(x,y)  
  3. * 初始条件为x=x[0]时,y=y[0].  
  4. * 输入: f--函数f(x,y)的指针  
  5. *       x--自变量离散值数组(其中x[0]为初始条件)  
  6. *       y--对应于自变量离散值的函数值数组(其中y[0]为初始条件)  
  7. *       h--计算步长  
  8. *       n--步数  
  9. * 输出: x为说求解的自变量离散值数组  
  10. *       y为所求解对应于自变量离散值的函数值数组  
  11. ************************************************************************/   
  12. double runge_kuta(double(*f)(double,double),double x[],   
  13.  double y[],double h,int n)   
  14. {   
  15.     int i;   
  16.     double xs,ys,xp,yp,dy;   
  17.     xs=x[0]+n*h;   
  18.     for(i=0;i<n;i++)   
  19.     {   
  20.         ys=y;   
  21.         dy=(*f)(x,y); //k1    
  22.         y[i+1]=y+h*dy/6;   
  23.         xp=x+h/2;   
  24.         yp=ys+h*dy/2;   
  25.         dy=(*f)(xp,yp); //k2    
  26.         y[i+1]+=h*dy/3;   
  27.         yp=ys+h*dy/2;   
  28.         dy=(*f)(xp,yp);  //k3    
  29.         y[i+1]+=h*dy/3;   
  30.         xp+=h/2;   
  31.         yp=ys+h*dy;   
  32.         dy=(*f)(xp,yp); //k4    
  33.         y[i+1]+=h*dy/6;   
  34.         x[i+1]=xp;   
  35.         if(x[i+1]>=xs)   
  36.         return (0);   
  37.     }   
  38.     return(0);   
  39. }   
/************************************************************************ 
* 用四阶(定步长)龙格--库塔法求解初值问题,其中一阶微分方程未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、改进的欧拉方法求解微分方程初值问题

感觉精度比较低
  1. /************************************************************************  
  2. * 用改进的欧拉方法求解初值问题,其中一阶微分方程未y'=f(x,y)  
  3. * 初始条件为x=x[0]时,y=y[0].  
  4. * 输入: f--函数f(x,y)的指针  
  5. *       x--自变量离散值数组(其中x[0]为初始条件)  
  6. *       y--对应于自变量离散值的函数值数组(其中y[0]为初始条件)  
  7. *       h--计算步长  
  8. *       n--步数  
  9. * 输出: x为说求解的自变量离散值数组  
  10. *       y为所求解对应于自变量离散值的函数值数组  
  11. ************************************************************************/   
  12. double proved_euler(double(*f)(double,double),double x[],   
  13. double y[],double h,int n)   
  14. {   
  15.     int i;   
  16.     double xs,ys,yp;   
  17.   
  18.     for(i=0;i<n;i++)   
  19.     {   
  20.         ys=y;   
  21.         xs=x;   
  22.         y[i+1]=y;   
  23.         yp=(*f)(xs,ys); //k1    
  24.         y[i+1]+=yp*h/2.0;   
  25.         ys+=h*yp;   
  26.         xs+=h;   
  27.         yp=(*f)(xs,ys); //k2    
  28.         y[i+1]+=yp*h/2.0;   
  29.         x[i+1]=xs;   
  30.     }   
  31.     return(0);   
  32. }   
/************************************************************************ 
* 用改进的欧拉方法求解初值问题,其中一阶微分方程未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。 中心差分(矩形)公式求导
  1. /************************************************************************  
  2. * 中心差分(矩形)公式计算函数f(x)在a点的导数值  
  3. * 输入: f--函数f(x)的指针  
  4. *       a--求导点  
  5. *       h--初始步长  
  6. *       eps--计算精度  
  7. *       max_it--最大循环次数  
  8. * 输出: 返回值为f(x)在a点的导数  
  9. ************************************************************************/   
  10. double central_difference(double (*f)(double),double a,   
  11.  double h,double eps,int max_it)   
  12. {   
  13.     double ff,gg;   
  14.     int k;   
  15.     ff=0.0;   
  16.     for(k=0;k<max_it;k++)   
  17.     {   
  18.         gg=((*f)(a+h)-(*f)(a-h))/(h+h);   
  19.         if(fabs(gg-ff)<eps)   
  20.         return(gg);   
  21.         h*=0.5;   
  22.         ff=gg;   
  23.     }   
  24.     if(k==max_it)   
  25.     {   
  26.         printf("未能达到精度要求,需增大迭代次数!");   
  27.         return(0);   
  28.     }   
  29.     return(gg);   
  30. }   
/************************************************************************ 
* 中心差分(矩形)公式计算函数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点法求积分

  1. /********************************************************************  
  2. * 用高斯10点法计算函数f(x)从a到b的积分值  
  3. * 输入: f--函数f(x)的指针  
  4. *       a--积分下限  
  5. *       b--积分上限  
  6. * 输出: 返回值为f(x)从a点到b点的积分值  
  7. *******************************************************************/   
  8. double gauss_legendre(double(*f)(double),double a,double b)   
  9. {   
  10.     const int n=10;   
  11.     const double z[10]={-0.9739065285,-0.8650633677,-0.6794095683,   
  12.         -0.4333953941,-0.1488743390,0.1488743390,   
  13.         0.4333953941,0.6794095683,0.8650633677,   
  14.         0.9739065285};   
  15.     const double w[10]={0.0666713443,0.1494513492,0.2190863625,   
  16.         0.2692667193,0.2955242247,0.2955242247,   
  17.         0.2692667193,0.2190863625,0.1494513492,   
  18.         0.0666713443};   
  19.     double y,gg;   
  20.     int i;   
  21.     gg=0.0;   
  22.     for(i=0;i<n;i++)   
  23.     {   
  24.         y=(z[i]*(b-a)+a+b)/2.0;   
  25.         gg+=w[i]*(*f)((double)y);   
  26.     }   
  27.     return((double)((gg*(b-a)/2.0)));   
  28. }   
/******************************************************************** 
* 用高斯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、龙贝格法求积分

  1. /********************************************************************  
  2. * 用龙贝格法计算函数f(x)从a到b的积分值  
  3. * 输入: f--函数f(x)的指针  
  4. *       a--积分下限  
  5. *       b--积分上限  
  6. *       eps--计算精度  
  7. *       max_it--最大迭代次数  
  8. * 输出: 返回值为f(x)从a点到b点的积分值  
  9. *******************************************************************/   
  10. double romberg(double(*f)(double),double a,double b,   
  11.                           double eps,int max_it)   
  12. {   
  13.     double *t,h;   
  14.     int i,m,k;   
  15.   
  16.     if(!(t=(double *)malloc(max_it*sizeof(double)+1)))   
  17.         return(ERROR_CODE);   
  18.     h=b-a;   
  19.     t[1]=h*((*f)(a)+(*f)(b))/2.0;   
  20.     printf("%18.10e\n",t[1]);   
  21.     for(k=2;k<max_it+1;k++)   
  22.     {   
  23.         double s,sm;   
  24.         h*=0.5;   
  25.         s=0.0;   
  26.         for(i=0;i<pow(2,k-2);i++)   
  27.             s+=(*f)(a+(2*i+1)*h);   
  28.         sm=t[1];   
  29.         t[1]=0.5*t[1]+h*s;   
  30.         for(m=2;m<k+1;m++)   
  31.         {   
  32.             s=t[m];   
  33.             t[m]=t[m-1]+(t[m-1]-sm)/(pow(4,m-1)-1);   
  34.             if(m<k)   
  35.                 sm=s;   
  36.         }   
  37.         for(m=1;m<k+1;m++)   
  38.             printf("%18.10e",t[m]);   
  39.         printf("\n");   
  40.         if(fabs(t[k]-sm)<eps)   
  41.         {   
  42.             sm=t[k];   
  43.             free(t);   
  44.             return(sm);   
  45.         }   
  46.     }   
  47.     return(ERROR_CODE);   
  48. }   
/******************************************************************** 
* 用龙贝格法计算函数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、复合辛普生法求积分

  1. #include "stdio.h"    
  2.   
  3. double composite_simpson(double(*f)(double),double a,double b,int n);   
  4. double myfun(double x);   
  5.   
  6. void main()   
  7. {   
  8.     double(*fun)(double);   
  9.     double a,b,S4;   
  10.     int n;   
  11.   
  12.     a=0;   
  13.     b=1;   
  14.     n=4;   
  15.     fun=myfun;   
  16.     S4=composite_simpson(fun,a,b,n);   
  17.     printf("\n积分值为:%f\n",S4);   
  18. }   
  19.   
  20. /********************************************************************  
  21. * 用复合辛普生法计算函数f(x)从a到b的积分值  
  22. * 输入: f--函数f(x)的指针  
  23. *       a--积分下限  
  24. *       b--积分上限  
  25. *       n--分段数  
  26. * 输出: 返回值为f(x)从a点到b点的积分值  
  27. *******************************************************************/   
  28. double composite_simpson(double(*f)(double),double a,double b,int n)   
  29. {   
  30.     double s,h,x;   
  31.     int i;   
  32.     printf("x\t\tf(x)\t\ts\n");   
  33.     s=(*f)(a)-(*f)(b);   
  34.     h=(b-a)/(2*n);   
  35.     x=a;   
  36.     for(i=1;i<2*n;i+=2)   
  37.     {   
  38.         x+=h;   
  39.         s+=4*(*f)(x);   
  40.         printf("%f\t%f\t%f\n",x,(*f)(x),s*h/3);   
  41.         x+=h;   
  42.         s+=2*(*f)(x);   
  43.         printf("%f\t%f\t%f\n",x,(*f)(x),s*h/3);   
  44.     }   
  45.     return(s*h/3);   
  46. }   
  47.   
  48.   
  49. double myfun(double x)   
  50. {   
  51.     double y;   
  52.   
  53.     y=4.0/(1+x*x);   
  54.   
  55.     return(y);   
  56. }   
#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、最小二乘法拟合

  1. /***************************************************************  
  2. * 本算法用最小二乘法依据指定的M个基函数及N个已知数据进行曲线拟和  
  3. * 输入: m--已知数据点的个数M  
  4. *       f--M维基函数向量  
  5. * n--已知数据点的个数N-1  
  6. *       x--已知数据点第一坐标的N维列向量  
  7. *       y--已知数据点第二坐标的N维列向量  
  8. *       a--无用  
  9. * 输出: 函数返回值为曲线拟和的均方误差  
  10. *       a为用基函数进行曲线拟和的系数,  
  11. *       即a[0]f[0]+a[1]f[1]+...+a[M]f[M].  
  12. ****************************************************************/   
  13. double mini_product(int m,double(*f[M])(double),int n,double x[N],   
  14. double y[N],double a[M])   
  15. {   
  16.     double e,ff,b[M][M],c[M][1];   
  17.     int i,j,k;   
  18.   
  19.     for(j=0;j<m;j++)       /*计算最小均方逼近矩阵及常向量*/   
  20.     {   
  21.         for(k=0;k<m;k++)   
  22.         {   
  23.             b[j][k]=0.0;   
  24.             for(i=0;i<n;i++)   
  25.                 b[j][k]+=(*f[j])(x)*(*f[k])(x);   
  26.         }   
  27.         c[j][0]=0.0;   
  28.         for(i=0;i<n;i++)   
  29.             c[j][0]+=(*f[j])(x)*y;   
  30.     }   
  31.     gaussian_elimination(m,b,1,c);   /*求拟和系数*/   
  32.     for(i=0;i<m;i++)   
  33.     a=c[0];   
  34.     e=0.0;   
  35.     for(i=0;i<n;i++)   /*计算均方误差*/   
  36.     {   
  37.         ff=0.0;   
  38.         for(j=0;j<m;j++)   
  39.             ff+=a[j]*(*f[j])(x);   
  40.         e+=(y-ff)*(y-ff);   
  41.     }   
  42.     return(e);   
  43. }   
  44.   
  45. /*************************************************************************  
  46. * 高斯列主元素消去法求解矩阵方程AX=B,其中A是N*N的矩阵,B是N*M矩阵  
  47. * 输入: n----方阵A的行数  
  48. *       a----矩阵A  
  49. *       m----矩阵B的列数  
  50. *       b----矩阵B  
  51. * 输出: det----矩阵A的行列式值  
  52. *       a----A消元后的上三角矩阵  
  53. *       b----矩阵方程的解X  
  54. **************************************************************************/   
  55. double gaussian_elimination(int n,double a[M][M],int m,double b[M][1])   
  56. {   
  57.     int i,j,k,mk;   
  58.     double det,mm,f;   
  59.     det = 1.0;   
  60.     for(k = 0;k<n-1;k++)  /*选主元并消元*/   
  61.     {   
  62.         mm=a[k][k];   
  63.         mk = k;   
  64.         for(i=k+1;i<n;i++)   /*选择第K列主元素*/   
  65.         {   
  66.             if(fabs(mm)<fabs(a<EM>[k]))   
  67.             {   
  68.                 mm = a<EM>[k];   
  69.                 mk = i;   
  70.             }   
  71.         }   
  72.         if(fabs(mm)<EPS)   
  73.             return(0);   
  74.         if(mk!=k) /* 将第K列主元素换行到对角线上*/   
  75.         {   
  76.             for(j=k;j<n;j++)   
  77.             {   
  78.                 f = a[k][j];   
  79.                 a[k][j]=a[mk][j];   
  80.                 a[mk][j]=f;   
  81.             }   
  82.             for(j=0;j<m;j++)   
  83.             {   
  84.                 f = b[k][j];   
  85.                 b[k][j]=b[mk][j];   
  86.                 b[mk][j]=f;   
  87.             }   
  88.             det = -det;   
  89.         }   
  90.         for(i=k+1;i<n;i++) /*将第K列对角线以下消元为零*/   
  91.         {   
  92.             mm = a<EM>[k]/a[k][k];   
  93.             a<EM>[k]=0.0;   
  94.             for(j=k+1;j<n;j++)   
  95.                 a<EM>[j]=a<EM>[j]-mm*a[k][j];   
  96.             for(j=0;j<m;j++)   
  97.                 b<EM>[j]=b<EM>[j]-mm*b[k][j];   
  98.         }   
  99.         det = det*a[k][k];   
  100.     }   
  101.     if(fabs(a[k][k])<EPS)   
  102.         return 0;   
  103.     det=det*a[k][k];   
  104.     for(i=0;i<m;i++) /*回代求解*/   
  105.     {   
  106.         b[n-1]<EM>=b[n-1]<EM>/a[n-1][n-1];   
  107.         for(j=n-2;j>=0;j--)   
  108.         {   
  109.             for(k=j+1;k<n;k++)   
  110.             b[j]<EM>=b[j]<EM>-a[j][k]*b[k]<EM>;   
  111.             b[j]<EM>=b[j]<EM>/a[j][j];   
  112.         }   
  113.     }   
  114.     return(det);   
  115. } </EM></EM></EM></EM></EM></EM></EM></EM></EM></EM></EM></EM></EM></EM></EM>  
/*************************************************************** 
* 本算法用最小二乘法依据指定的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.埃特金插值法

  1. /******************************************************  
  2. * 用埃特金插值法依据N个已知数据点计算函数值  
  3. * 输入: n--已知数据点的个数N-1  
  4. *       x--已知数据点第一坐标的N维列向量  
  5. * y--已知数据点第二坐标的N维列向量  
  6. * xx-插值点第一坐标  
  7. *       eps--求解精度  
  8. * 输出: 函数返回值所求插值点的第二坐标  
  9. ******************************************************/   
  10. double aitken(int n,double x[N],double y[N],double xx,double eps)   
  11. {   
  12.     double d[N];   
  13.     int i,j;   
  14.   
  15.     for(i=0;i<=n;i++)   
  16.         d=y;   
  17.     for(i=0;i<=n;i++)   
  18.     {   
  19.         for(j=0;j<i;j++)   
  20.             d=(d*(x[j]-xx)-d[j]*(x-xx))/(x[j]-x);   
  21.         if(d-d[i-1]<eps)   
  22.         return(d);   
  23.     }   
  24. }   
/****************************************************** 
* 用埃特金插值法依据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、牛顿插值法

  1. /******************************************************  
  2. * 用牛顿插值法依据N个已知数据点即使函数值  
  3. * 输入: n--已知数据点的个数N-1  
  4. *       x--已知数据点第一坐标的N维列向量  
  5. * y--已知数据点第二坐标的N维列向量  
  6. * xx-插值点第一坐标  
  7. * 输出: 函数返回值所求插值点的第二坐标  
  8. ******************************************************/   
  9. double newton(int n,double x[N],double y[N],double xx)   
  10. {   
  11.     double d[N],b;   
  12.     int i,j;   
  13.   
  14.     for(i=0;i<=n;i++)   
  15.         d=y;   
  16.     for(i=n-1;i>=0;i--) /*求差商*/   
  17.         for(j=i+1;j<=n;j++)   
  18.         {   
  19.             if(fabs(x-x[j])<EPS)   
  20.                 return 0;   
  21.             d[j]=(d[j-1]-d[j])/(x-x[j]);   
  22.         }   
  23.     b=d[n];   
  24.     for(i=n-1;i>=0;i--)   
  25.     b=d+(xx-x)*b;   
  26.     return b;   
  27. }   
/****************************************************** 
* 用牛顿插值法依据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、拉格朗日插值法

  1. /******************************************************  
  2. * 用拉格朗日插值法依据N个已知数据点即使函数值  
  3. * 输入: n--已知数据点的个数N-1  
  4. *       x--已知数据点第一坐标的N维列向量  
  5. * y--已知数据点第二坐标的N维列向量  
  6. * xx-插值点第一坐标  
  7. * 输出: 函数返回值所求插值点的第二坐标  
  8. ******************************************************/   
  9. double lagrange(int n,double x[N],double y[N],double xx)   
  10. {   
  11.     double p,yy;   
  12.     int i,j;   
  13.     yy = 0.0;   
  14.     for(i=0;i<=n;i++)   
  15.     {   
  16.         p=1.0;   
  17.         for(j=0;j<=n;j++)   
  18.             if(i!=j)   
  19.             {   
  20.                 if(fabs(x-x[j])<EPS)   
  21.                     return 0;   
  22.                 p=p*(xx-x[j])/(x-x[j]);   
  23.             }   
  24.         yy=yy+p*y;   
  25.     }   
  26.     return(yy);   
  27. }   
/****************************************************** 
* 用拉格朗日插值法依据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

  1. /*************************************************************************  
  2. * 逆矩阵法求解矩阵方程AX=B,其中A是N*N的矩阵,B是N*N矩阵  
  3. * 输入: n----方阵A的行数  
  4. *       a----矩阵A  
  5. *       m----矩阵B的列数  
  6. *       b----矩阵B  
  7. * 输出: det----矩阵A的行列式值  
  8. *       a----A的逆矩阵  
  9. *       b----矩阵方程的解X  
  10. **************************************************************************/   
  11. double gaussian_jodan_solve(int n,double a[N][N],int m,double b[N][M])   
  12. {   
  13.     double det,f[N];   
  14. int i,j,k;   
  15.       
  16.     det = gaussian_jodan(n,a);   
  17.     if(det==0) return (0);   
  18.     for(k=0;k<m;k++) /*做矩阵乘法AB*/   
  19.     {   
  20.         for(i=0;i<n;i++)   
  21.         {   
  22.             f=0.0;   
  23.             for(j=0;j<n;j++)   
  24.                 f=f+a[j]*b[j][k];   
  25.         }   
  26.         for(i=0;i<n;i++)   
  27.             b[k]=f;   
  28.     }   
  29.     return(det);   
  30. }   
/************************************************************************* 
* 逆矩阵法求解矩阵方程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); 
} 


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

  1. /*************************************************************************  
  2. * 高斯--约当列主元素法求矩阵方程A的逆矩阵,其中A是N*N的矩阵  
  3. * 输入: n----方阵A的行数  
  4. *       a----矩阵A  
  5. * 输出: det--A的行列式的值  
  6. *       a----A的逆矩阵  
  7. **************************************************************************/   
  8. double gaussian_jodan(int n,double a[N][N])   
  9. {   
  10.     int i,j,k,mk;   
  11.     int p[N];  /*记录主行元素在原矩阵中的位置*/   
  12.     double det,m,f;   
  13.   
  14.     det = 1.0;   
  15.     for(k=0;k<n;k++)   
  16.     {   
  17.         m=a[k][k];  /*选第K列主元素*/   
  18.         mk=k;   
  19.         for(i=k+1;i<n;i++)   
  20.             if(fabs(m)<fabs(a[k]))   
  21.             {   
  22.                 m=a[k];   
  23.                 mk=i;   
  24.             }   
  25.         if(fabs(m)<EPS) return(0);   
  26.         if(mk!=k)   
  27.         {   
  28.             for(j=0;j<n;j++) /*将第K列主元素换行到主对角线上*/   
  29.             {   
  30.                 f=a[k][j];   
  31.                 a[k][j]=a[mk][j];   
  32.                 a[mk][j]=f;   
  33.             }   
  34.             p[k]=mk;   
  35.             det = -det;   
  36.         }   
  37.         else   
  38.             p[k]=k;   
  39.         det=det*m;   
  40.         for(j=0;j<n;j++)  /*计算主行元素*/   
  41.             if(j!=k)   
  42.                 a[k][j]=a[k][j]/a[k][k];   
  43.         a[k][k]=1.0/a[k][k];   
  44.         for(i=0;i<n;i++) /*消元*/   
  45.         {   
  46.             if(i!=k)   
  47.             {   
  48.                 for(j=0;j<n;j++)   
  49.                     if(j!=k)   
  50.                 a[j]=a[j]-a[k]*a[k][j];   
  51.                 a[k]=-a[k]*a[k][k];   
  52.             }   
  53.         }   
  54.     }   
  55.     for(k=n-2;k>=0;k--) /*按主行在原矩阵中的位置交换列*/   
  56.     {   
  57.         if(p[k]!=k)   
  58.         for(i=0;i<n;i++)   
  59.         {   
  60.             f=a[k];   
  61.             a[k]=a[p[k]];   
  62.             a[p[k]]=f;   
  63.         }   
  64.     }   
  65.     return(det);   
  66. }   
/************************************************************************* 
* 高斯--约当列主元素法求矩阵方程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、高斯列主元素消去法求解矩阵方程

  1. /*************************************************************************  
  2. * 高斯列主元素消去法求解矩阵方程AX=B,其中A是N*N的矩阵,B是N*M矩阵  
  3. * 输入: n----方阵A的行数  
  4. *       a----矩阵A  
  5. *       m----矩阵B的列数  
  6. *       b----矩阵B  
  7. * 输出: det----矩阵A的行列式值  
  8. *       a----A消元后的上三角矩阵  
  9. *       b----矩阵方程的解X  
  10. **************************************************************************/   
  11. double gaussian_elimination(int n,double a[N][N],int m,double b[N][M])   
  12. {   
  13.     int i,j,k,mk;   
  14.     double det,mm,f;   
  15.     det = 1.0;   
  16.     for(k = 0;k<n-1;k++)  /*选主元并消元*/   
  17.     {   
  18.         mm=a[k][k];   
  19.         mk = k;   
  20.         for(i=k+1;i<n;i++)   /*选择第K列主元素*/   
  21.         {   
  22.             if(fabs(mm)<fabs(a[k]))   
  23.             {   
  24.                 mm = a[k];   
  25.                 mk = i;   
  26.             }   
  27.         }   
  28.         if(fabs(mm)<EPS)   
  29.             return(0);   
  30.         if(mk!=k) /* 将第K列主元素换行到对角线上*/   
  31.         {   
  32.             for(j=k;j<n;j++)   
  33.             {   
  34.                 f = a[k][j];   
  35.                 a[k][j]=a[mk][j];   
  36.                 a[mk][j]=f;   
  37.             }   
  38.             for(j=0;j<m;j++)   
  39.             {   
  40.                 f = b[k][j];   
  41.                 b[k][j]=b[mk][j];   
  42.                 b[mk][j]=f;   
  43.             }   
  44.             det = -det;   
  45.         }   
  46.         for(i=k+1;i<n;i++) /*将第K列对角线以下消元为零*/   
  47.         {   
  48.             mm = a[k]/a[k][k];   
  49.             a[k]=0.0;   
  50.             for(j=k+1;j<n;j++)   
  51.                 a[j]=a[j]-mm*a[k][j];   
  52.             for(j=0;j<m;j++)   
  53.                 b[j]=b[j]-mm*b[k][j];   
  54.         }   
  55.         det = det*a[k][k];   
  56.     }   
  57.     if(fabs(a[k][k])<EPS)   
  58.         return 0;   
  59.     det=det*a[k][k];   
  60.     for(i=0;i<m;i++) /*回代求解*/   
  61.     {   
  62.         b[n-1]=b[n-1]/a[n-1][n-1];   
  63.         for(j=n-2;j>=0;j--)   
  64.         {   
  65.             for(k=j+1;k<n;k++)   
  66.                 b[j]=b[j]-a[j][k]*b[k];   
  67.             b[j]=b[j]/a[j][j];   
  68.         }   
  69.     }   
  70.     return(det);   
  71. }   
/************************************************************************* 
* 高斯列主元素消去法求解矩阵方程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乘以一个象素的面积就是所求的面积了。  

http://www.cnblogs.com/haohello/archive/2007/04/26/727744.html
  • 2
    点赞
  • 10
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值