数值分析多种算法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、四阶(定步长)龙格--库塔法求解微分初值问题
精度比欧拉方法高
但是感觉依然不理想
|
/************************************************************************ |
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++) |
7。中心差分(矩形）公式求导
|
/************************************************************************ |
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))) |
10、复合辛普生法求积分
|
#include "stdio.h" double composite_simpson(double(*f)(double),double a,double b,int n); void main() a=0; } /******************************************************************** } double myfun(double x) y=4.0/(1+x*x); return(y); } |
11、最小二乘法拟合
|
/*************************************************************** for(j=0;j<m;j++) /*计算最小均方逼近矩阵及常向量*/ |
|
/************************************************************************* * 高斯列主元素消去法求解矩阵方程AX=B,其中A是N*N的矩阵,B是N*M矩阵 * 输入: n----方阵A的行数 * a----矩阵A * m----矩阵B的列数 * b----矩阵B * 输出: det----矩阵A的行列式值 * a----A消元后的上三角矩阵 * b----矩阵方程的解X **************************************************************************/ double gaussian_elimination(int n,double a[M][M],int m,double b[M][1]) { int i,j,k,mk; double det,mm,f; det = 1.0; for(k = 0;k<n-1;k++) /*选主元并消元*/ { mm=a[k][k]; mk = k; for(i=k+1;i<n;i++) /*选择第K列主元素*/ { if(fabs(mm)<fabs(a[k])) { mm = a[k]; mk = i; } } if(fabs(mm)<EPS) return(0); if(mk!=k) /* 将第K列主元素换行到对角线上*/ { for(j=k;j<n;j++) { f = a[k][j]; a[k][j]=a[mk][j]; a[mk][j]=f; } for(j=0;j<m;j++) { f = b[k][j]; b[k][j]=b[mk][j]; b[mk][j]=f; } det = -det; } for(i=k+1;i<n;i++) /*将第K列对角线以下消元为零*/ { mm = a[k]/a[k][k]; a[k]=0.0; for(j=k+1;j<n;j++) a[j]=a[j]-mm*a[k][j]; for(j=0;j<m;j++) b[j]=b[j]-mm*b[k][j]; } det = det*a[k][k]; } if(fabs(a[k][k])<EPS) return 0; det=det*a[k][k]; for(i=0;i<m;i++) /*回代求解*/ { b[n-1]=b[n-1]/a[n-1][n-1]; for(j=n-2;j>=0;j--) { for(k=j+1;k<n;k++) b[j]=b[j]-a[j][k]*b[k]; b[j]=b[j]/a[j][j]; } } return(det); } |
12.埃特金插值法
|
/****************************************************** * 用埃特金插值法依据N个已知数据点计算函数值 * 输入: n--已知数据点的个数N-1 * x--已知数据点第一坐标的N维列向量 * y--已知数据点第二坐标的N维列向量 * xx-插值点第一坐标 * eps--求解精度 * 输出: 函数返回值所求插值点的第二坐标 ******************************************************/ double aitken(int n,double x[N],double y[N],double xx,double eps) { double d[N]; int i,j; for(i=0;i<=n;i++) |
13、牛顿插值法
|
/****************************************************** for(i=0;i<=n;i++) |
14、 拉格朗日插值法
|
/****************************************************** * 用拉格朗日插值法依据N个已知数据点即使函数值 * 输入: n--已知数据点的个数N-1 * x--已知数据点第一坐标的N维列向量 * y--已知数据点第二坐标的N维列向量 * xx-插值点第一坐标 * 输出: 函数返回值所求插值点的第二坐标 ******************************************************/ double lagrange(int n,double x[N],double y[N],double xx) { double p,yy; int i,j; yy = 0.0; for(i=0;i<=n;i++) { p=1.0; for(j=0;j<=n;j++) if(i!=j) { if(fabs(x-x[j])<EPS) return 0; p=p*(xx-x[j])/(x-x[j]); } yy=yy+p*y; } return(yy); } |
15、 逆矩阵法求解矩阵方程AX=B
|
/************************************************************************* det = gaussian_jodan(n,a); |
调用到的求逆矩阵的子函数
|
/************************************************************************* * 高斯--约当列主元素法求矩阵方程A的逆矩阵，其中A是N*N的矩阵 * 输入: n----方阵A的行数 * a----矩阵A * 输出: det--A的行列式的值 * a----A的逆矩阵 **************************************************************************/ double gaussian_jodan(int n,double a[N][N]) { int i,j,k,mk; int p[N]; /*记录主行元素在原矩阵中的位置*/ double det,m,f; det = 1.0; |
16、高斯列主元素消去法求解矩阵方程
|
/************************************************************************* |
17、任意多边形的面积计算（包括凹多边形的）
任意多边形的面积计算（包括凹多边形的，以及画多边形时线相交了的判别），最好提供相关资料，越详细越好，谢谢
---------------------------------------------------------------
我们都知道已知A（x1,y1）B（x2,y2）C（x3,y3）三点的面积公式为
¦x1 x2 x3 ¦
S(A,B,C) = ¦y1 y2 y3 ¦ * 0.5 (当三点为逆时针时为正，顺时针则为负的)
¦1 1 1 ¦
对多边形A1A2A3、、、An（顺或逆时针都可以），设平面上有任意的一点P，则有：
S（A1，A2，A3，、、、，An）
= abs(S（P，A1，A2） + S（P，A2，A3）+、、、+S（P，An，A1）)
P是可以取任意的一点，用（0，0）就可以了。
这种方法对凸和凹多边形都适用。
还有一个方法：
任意一个简单多边形，当它的各个顶点位于网格的结点上时，它的面积数S=b/2+c+1
其中：b代表该多边形边界上的网络结点数目
c代表该多边形内的网络结点数目
所以把整个图形以象素为单位可以把整个图形分成若干个部分，计算该图形边界上的点b和内部的点c就得到面积数S了，然后把S乘以一个象素的面积就是所求的面积了。
非常感谢“山城棒棒儿的MATLAB&FPGA世界 ”，从中获益匪浅，有时间也会整理一些相关的数值分析的代码，共享给急切需要，没有目标的网友同行们！希望能在您能获多获少都会有所收获，那将是我最幸福的事！
http://www.cnblogs.com/haohello/archive/2007/04/26/727744.html
数值分析多种算法C语言代码相关推荐
- 数值分析多种算法C语言代码－推荐
1.离散傅立叶变换与反变换 /************************************************************************ * 离散傅立叶变换与反变 ...
- 匈牙利算法c语言代码,漫谈匈牙利算法
匈牙利算法是由匈牙利数学家Edmonds于1965年提出,因而得名.匈牙利算法是基于Hall定理中充分性证明的思想,它是部图匹配最常见的算法,该算法的核心就是寻找增广路径,它是一种用增广路径求二分图最 ...
- 统计学习导论之R语言应用（四）：分类算法R语言代码实战
统计学习导论之R语言应用(ISLR) 参考资料: The Elements of Statistical Learning An Introduction to Statistical Learnin ...
- 相关数值分析多种算法代码
整理一些相关的数值分析的代码,共享给急切需要同行们!希望能在您能获多获少都会有所收获.>_<呵呵. 离散傅立叶变换与反变换 //****************************** ...
- 均匀不变lbp算法c语言代码,LBP原理介绍以及算法实现
有些读者可能会觉得奇怪,上次推送怎么就突然说起了双线性插值而不是继续介绍经典人脸识别的算法,其实是因为在学习圆形LBP算子的时候发现需要用到双线性插值于是顺带介绍了一下.(因为个人原因没经常看公众号的 ...
- 最小生成树普里姆算法c语言代码,普里姆算法生成最小生成树-C语言描述.doc
PAGE JIN JINGCHU UNIVERSITY OF TECHNOLOGY <数据结构(C语言描述)> 课程设计 学 院 计算机工程学院 班 级 12级软件技术1班 学 号 201 ...
- GSM A5/1算法C语言代码实现和分析
介绍 全球超过200个国家和地区超过10亿人正在使用GSM电话.对中国用户来说,GSM就是移动和联通的2g模式. 在1982年A5首次提出时,人们认为A5 / 1密钥长度要128位,但最终确定的结果是 ...
- pso算法c++语言代码,一C++PSO（PSO）算法
收集和变化PSO算法,它可用于参考实施: #include #include #include #include #include #define rand_01 ((float)rand() / ( ...
- OPT和LRU页面置换算法C语言代码,页面置换算法模拟——OPT、FIFO和LRU算法.doc
实用标准文案 精彩文档 操作系统实验报告 页面置换算法模拟 --OFT.FIFO和LRU算法 班级:2013级软件工程1班 学号:X X X 姓名:萧氏一郎 数据结构说明: Memery[10]物理块 ...
最新文章
- 全面访问JavaScript的最佳资源
- 工具安装===Sublime Text-安装
- JavaScript实现冒泡排序 可视化
- python爬虫天气预报_Python爬虫实例扒取2345天气预报
- Android--多个Activity共享Socket--单例模式
- 使用BUCK进行iOS项目打包
- 每个程序员都该学习的5种开发语言，不可错过！
- java8如何兼容java7_尽管使用Java 8功能,项目如何支持Java 7
- 05:年龄与疾病【一维数组】
- 编译指令#pragma详解
- sql: sql developer tunnel转接
- 3套看漫画学python视频教程
- 河北省人民检察院利用深信服桌面云办公，实现智慧检务
- 数学基础知识总结 —— 2. 常用积分公式
- 而立之年才感悟到的一些箴言:
- Java file.encoding
- Gephi简单导入数据
- linux 字体显示更清晰,Fedora下使中文字体显示变得更清晰
- python智力问答游戏代码,python实现智力问答测试小程序
- 普渡大学计算机图形,普渡大学计算机图形学技术研究生语言及申请要求-费用-课程设置...