数值分析C代码.docx

1、数值分析C代码目录1、离散傅立叶变换与反变换 12、四阶亚当姆斯预估计求解初值问题 23、几种常见随机数的产生 34、指数平滑法预测数据 85、四阶(定步长)龙格-库塔法求解微分初值问题 96、改进的欧拉方法求解微分方程初值问题 107、中心差分(矩形）公式求导 128、高斯10点法求积分 139、龙贝格法求积分 1410、复合辛普生法求积分 1511、最小二乘法拟合 1712.埃特金插值法 2113、牛顿插值法 2214、拉格朗日插值法 2315、逆矩阵法求解矩阵方程AX=B 2416、高斯列主元素消去法求解矩阵方程 2817、任意多边形的面积计算（包括凹多边形的） 301、离散傅立叶变换与

2、反变换/* 离散傅立叶变换与反变换* 输入: 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;kn;k+)w=k*q;ak=bk=0.0;for(i=0;in;i+)d=i*w;c=cos(d);s=sin(d)*sign;

3、ak+=c*x+s*y;bk+=c*y-s*x;if(sign=-1)c=1.0/n;for(k=0;kn;k+)ak=c*ak;bk=c*bk;2、四阶亚当姆斯预估计求解初值问题/* 用四阶亚当姆斯预估计求解初值问题，其中一阶微分方程未y=f(x,y)* 初始条件为x=x0时,yy0.* 输入: f-函数f(x,y)的指针* x-自变量离散值数组(其中x0为初始条件)* y-对应于自变量离散值的函数值数组(其中y0为初始条件)* h-计算步长* n-步数* 输出: x为说求解的自变量离散值数组* y为所求解对应于自变量离散值的函数值数组*/double adams(double(*f)(do

4、uble,double),double x,double y,double h,int n)double dy4,c,p,c1,p1,m;int i,j;runge_kuta(f,x,y,h,3);for(i=0;i4;i+)dy=(*f)(x,y);c=0.0; p=0.0;for(i=4;in+1;i+)x=xi-1+h;p1=yi-1+h*(55*dy3-59*dy2+37*dy1-9*dy0)/24;m=p1+251*(c-p)/270;c1=yi-1+h*(9*(*f)(x,m)+19*dy3-5*dy2+dy1)/24;y=c1-19*(c1-p1)/270;c=c1; p=p1;

5、for(j=0;j3;j+)dyj=dyj+1;dy3=(*f)(x,y);return(0);3、几种常见随机数的产生#include stdlib.h#include stdio.h#include math.hdouble 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);do

6、uble 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;i10;i+)for(j=0;j5

7、;j+)x=poisson(a,&s);mean+=x;printf(%-13.7f,x);printf(n);mean/=50;printf(平均值为:%-13.7fn,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)/

8、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;i12;i+)x+=uniform(0.0,1.0,seed);x=x-6.0;y=mean+x*sigma;return(y);/* 指数分布* 输入: beta-指数分布均值* seed-种子*/

9、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);elsex=beta*log(u2);retu

10、rn(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

11、 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=a);x=i-1;return(x);4、指数平滑法预测数据/* 本算法用指数平滑法预测数据* 输入: k-平滑周期* n-原始数据个数* m-预测步数* alfa-加权系数* x-指向原始数据数组指针* 输出: s1-返回值为指向一次平滑结果数组指针* s2-返回值为指向二次指数平滑结果数组指针* s3-返

12、回值为指向三次指数平滑结果数组指针* xx-返回值为指向预测结果数组指针*/void phyc(int k,int n,int m,double alfa,double xN_MAX,double s1N_MAX,double s2N_MAX,double s3N_MAX,double xxN_MAX)double a,b,c,beta;int i;s1k-1=0;for(i=0;ik;k+)s1k-1+=x;s1k-1/=k;for(i=k;i=n;i+)s1=alfa*x+(1-alfa)*s1i-1;s22*k-2=0;for(i=k-1;i2*k-1;i+)s22*k-2+=s1;s2

13、2*k-2/=k;for(i=2*k-1;i=n;i+)s2=alfa*s1+(1-alfa)*s2i-1;s33*k-3=0;for(i=2*k-2;i3*k-2;i+)s33*k-3+=s2;s33*k-3/=k;for(i=3*k-2;i=n;i+)s3=alfa*s2+(1-alfa)*s3i-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

14、+b*m+c*m*m;5、四阶(定步长)龙格-库塔法求解微分初值问题精度比欧拉方法高但是感觉依然不理想/* 用四阶(定步长)龙格-库塔法求解初值问题，其中一阶微分方程未y=f(x,y)* 初始条件为x=x0时,yy0.* 输入: f-函数f(x,y)的指针* x-自变量离散值数组(其中x0为初始条件)* y-对应于自变量离散值的函数值数组(其中y0为初始条件)* h-计算步长* n-步数* 输出: x为说求解的自变量离散值数组* y为所求解对应于自变量离散值的函数值数组*/double runge_kuta(double(*f)(double,double),double x,double y

15、,double h,int n)int i;double xs,ys,xp,yp,dy;xs=x0+n*h;for(i=0;i=xs)return (0);return(0);6、改进的欧拉方法求解微分方程初值问题感觉精度比较低/* 用改进的欧拉方法求解初值问题，其中一阶微分方程未y=f(x,y)* 初始条件为x=x0时,yy0.* 输入: f-函数f(x,y)的指针* x-自变量离散值数组(其中x0为初始条件)* y-对应于自变量离散值的函数值数组(其中y0为初始条件)* h-计算步长* n-步数* 输出: x为说求解的自变量离散值数组* y为所求解对应于自变量离散值的函数值数组*/doub

16、le proved_euler(double(*f)(double,double),double x,double y,double h,int n)int i;double xs,ys,yp;for(i=0;in;i+)ys=y;xs=x;yi+1=y;yp=(*f)(xs,ys); /k1yi+1+=yp*h/2.0;ys+=h*yp;xs+=h;yp=(*f)(xs,ys); /k2yi+1+=yp*h/2.0;xi+1=xs;return(0);7、中心差分(矩形）公式求导/* 中心差分(矩形）公式计算函数f(x)在a点的导数值* 输入: f-函数f(x)的指针* a-求导点* h-初

17、始步长* 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;kmax_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