数学算法程序.docx
- 文档编号:755744
- 上传时间:2023-04-30
- 格式:DOCX
- 页数:25
- 大小:22.43KB
数学算法程序.docx
《数学算法程序.docx》由会员分享,可在线阅读,更多相关《数学算法程序.docx(25页珍藏版)》请在冰点文库上搜索。
数学算法程序
1、离散傅立叶变换与反变换
/************************************************************************
*离散傅立叶变换与反变换
*输入:
x--要变换的数据的实部
*y--要变换的数据的虚部
* a--变换结果的实部
* b--变换结果的虚部
* n--数据长度
* sign--sign=1时,计算离散傅立叶正变换;sign=-1时;计算离散傅立叶反变换
************************************************************************/
voiddft(doublex[],doubley[],doublea[],doubleb[],intn,intsign)
{
inti,k;
doublec,d,q,w,s;
q=6.28318530718/n;
for(k=0;k { w=k*q; a[k]=b[k]=0.0; for(i=0;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 { 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为所求解对应于自变量离散值的函数值数组 ************************************************************************/ doubleadams(double(*f)(double,double),doublex[], doubley[],doubleh,intn) { doubledy[4],c,p,c1,p1,m; inti,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 { 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" doubleuniform(doublea,doubleb,longint*seed); doublegauss(doublemean,doublesigma,longint*seed); doubleexponent(doublebeta,longint*seed); doublelaplace(doublebeta,longint*seed); doublerayleigh(doublesigma,longint*seed); doubleweibull(doublea,doubleb,longint*seed); intbn(doublep,longint*seed); intbin(intn,doublep,longint*seed); intpoisson(doublelambda,longint*seed); voidmain() { doublea,b,x,mean; inti,j; longints; 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为随机数的种子 ********************************************************************/ doubleuniform(doublea,doubleb,longint*seed) { doublet; *seed=2045*(*seed)+1; *seed=*seed-(*seed/1048576)*1048576; t=(*seed)/1048576.0; t=a+(b-a)*t; return(t); } /******************************************************************* *正态分布 *输入: mean--双精度实型变量,正态分布的均值 * sigma--双精度实型变量,正态分布的均方差 * seed--长整型指针变量,*seed为随机数的种子 ********************************************************************/ doublegauss(doublemean,doublesigma,longint*seed) { inti; doublex,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--种子 *******************************************************************/ doubleexponent(doublebeta,longint*seed) { doubleu,x; u=uniform(0.0,1.0,seed); x=-beta*log(u); return(x); } /******************************************************************* *拉普拉斯随机分布 *beta--拉普拉斯分布的参数 **seed--随机数种子 *******************************************************************/ doublelaplace(doublebeta,longint*seed) { doubleu1,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); } /******************************************************************** *瑞利分布 * ********************************************************************/ doublerayleigh(doublesigma,longint*seed) { doubleu,x; u=uniform(0.,1.,seed); x=-2.0*log(u); x=sigma*sqrt(x); return(x); } /************************************************************************/ /*韦伯分布 */ /************************************************************************/ doubleweibull(doublea,doubleb,longint*seed) { doubleu,x; u=uniform(0.0,1.0,seed); u=-log(u); x=b*pow(u,1.0/a); return(x); } /************************************************************************/ /*贝努利分布 */ /************************************************************************/ intbn(doublep,longint*seed) { intx; doubleu; u=uniform(0.0,1.0,seed); x=(u<=p)? 1: 0; return(x); } /************************************************************************/ /*二项式分布 */ /************************************************************************/ intbin(intn,doublep,longint*seed) { inti,x; for(x=0,i=0;i x+=bn(p,seed); return(x); } /************************************************************************/ /*泊松分布 */ /************************************************************************/ intpoisson(doublelambda,longint*seed) { inti,x; doublea,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--返回值为指向预测结果数组指针 ************************************************************************/ voidphyc(intk,intn,intm,doublealfa,doublex[N_MAX], doubles1[N_MAX],doubles2[N_MAX],doubles3[N_MAX],doublexx[N_MAX]) { doublea,b,c,beta; inti; s1[k-1]=0; for(i=0;i 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为所求解对应于自变量离散值的函数值数组 ************************************************************************/ doublerunge_kuta(double(*f)(double,double),doublex[], doubley[],doubleh,intn) { inti; doublexs,ys,xp,yp,dy; xs=x[0]+n*h; for(i=0;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为所求解对应于自变量离散值的函数值数组 ************************************************************************/ doubleproved_euler(double(*f)(double,double),doublex[], doubley[],doubleh,intn) { inti; doublexs,ys,yp; for(i=0;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点的导数 ************************************************************************/ doublecentral_difference(double(*f)(double),doublea, doubleh,doubleeps,intmax_it) { doubleff,gg; intk; ff=0.0; for(k=0;k { gg=((*f)(a+h)-(*f)(a-h))/(h+h); if(fabs(gg-ff) 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点的积分值 *******************************************************************/ doublegauss_legendre(double(*f)(double),doublea,doubleb) { constintn=10; constdoublez[10]={-0.9739065285,-0.8650633677,-0.6794095683, -0.4333953941,-0.1488743390,0.1488743390, 0.4333953941,0.6794095683,0.8650633677, 0.9739065285}; constdoublew[10]={0.0666713443,0.1494513492,0.2190863625, 0.2692667193,0.2955242247,0.2955242247, 0.2692667193,0.2190863625,0.1494513492,
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数学 算法 程序