c.txt

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

数值分析多种算法c语言代码
1、离散傅立叶变换与反变换

/************************************************************************ 
* 离散傅立叶变换与反变换 
* 输入: 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。四阶亚当姆斯预估计求解初值问题

/************************************************************************ 
* 用四阶亚当姆斯预估计求解初值问题，其中一阶微分方程未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、几种常见随机数的产生

#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、指数平滑法预测数据

/************************************************************************ 
* 本算法用指数平滑法预测数据 
* 输入: 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、四阶(定步长)龙格--库塔法求解微分初值问题

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




/************************************************************************ 
* 用四阶(定步长)龙格--库塔法求解初值问题，其中一阶微分方程未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、改进的欧拉方法求解微分方程初值问题

感觉精度比较低 

/************************************************************************ 
* 用改进的欧拉方法求解初值问题，其中一阶微分方程未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。中心差分(矩形）公式求导


 



/************************************************************************ 
* 中心差分(矩形）公式计算函数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点法求积分

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 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))); 
} 
[/code] 


9、龙贝格法求积分


 


/******************************************************************** 
* 用龙贝格法计算函数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) 
{