数学算法程序.docx

数学算法程序.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,