数值分析C代码.docx

数值分析C代码.docx

《数值分析C代码.docx》由会员分享，可在线阅读，更多相关《数值分析C代码.docx（26页珍藏版）》请在冰豆网上搜索。

数值分析C代码

目录

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

2、四阶亚当姆斯预估计求解初值问题2

3、几种常见随机数的产生3

4、指数平滑法预测数据8

5、四阶（定步长）龙格--库塔法求解微分初值问题9

6、改进的欧拉方法求解微分方程初值问题10

7、中心差分（矩形）公式求导12

8、高斯10点法求积分13

9、龙贝格法求积分14

10、复合辛普生法求积分15

11、最小二乘法拟合17

12.埃特金插值法21

13、牛顿插值法22

14、拉格朗日插值法23

15、逆矩阵法求解矩阵方程AX=B24

16、高斯列主元素消去法求解矩阵方程28

17、任意多边形的面积计算（包括凹多边形的）30

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） 

{ 

constint