三点高斯积分程序

原创 2007年09月17日 18:04:00
//高斯积分程序,高斯积分具有计算速度快,精度高,能计算反常积分等优点
#include <stdio.h>
#include <math.h>

#define PI 3.14159265
static const double gp3[3]={-0.77459666924148338, 0.0                , 0.77459666924148338};//高斯点
static const double gc3[3]={ 0.55555555555555556, 0.88888888888888889, 0.55555555555555556};//系数

//三点Guass积分
double Guass_Integral(double x1,double x2,double (*fun)(double),int n)
{
    int i;
    double dx,k,l,r=0.0;
    dx=(x2-x1)/n;
    k=dx/2;
    for(i=0;i<n;i++)
    {
        l=x1+(i+0.5)*dx;
        r+=gc3[0]*fun(k*gp3[0]+l)+gc3[1]*fun(k*gp3[1]+l)+gc3[2]*fun(k*gp3[2]+l);
    }
    return r*k;
}

double Guass_Integral(double u1,double u2,double v1,double v2,void (*fun)(double,double,double *,double *,double *),int m,int n)
{
    int i,j;
    double rx=0.0,ry=0.0,rz=0.0;
    double x00,x01,x02,x10,x11,x12,x20,x21,x22;
    double y00,y01,y02,y10,y11,y12,y20,y21,y22;
    double z00,z01,z02,z10,z11,z12,z20,z21,z22;
    double du,dv,ku,kv,lu,lv;
    du=(u2-u1)/m;
    dv=(v2-v1)/n;
    ku=du/2;
    kv=dv/2;
    for(i=0;i<m;i++)
    {
        for(j=0;j<n;j++)
        {
            lu=u1+(i+0.5)*du;
            lv=v1+(j+0.5)*dv;
            //为了计算速度快,展开循环
            fun(ku*gp3[0]+lu,kv*gp3[0]+lv,&x00,&y00,&z00);
            fun(ku*gp3[0]+lu,kv*gp3[1]+lv,&x01,&y01,&z01);
            fun(ku*gp3[0]+lu,kv*gp3[2]+lv,&x02,&y02,&z02);
            fun(ku*gp3[1]+lu,kv*gp3[0]+lv,&x10,&y10,&z10);
            fun(ku*gp3[1]+lu,kv*gp3[1]+lv,&x11,&y11,&z11);
            fun(ku*gp3[1]+lu,kv*gp3[2]+lv,&x12,&y12,&z12);
            fun(ku*gp3[2]+lu,kv*gp3[0]+lv,&x20,&y20,&z20);
            fun(ku*gp3[2]+lu,kv*gp3[1]+lv,&x21,&y21,&z21);
            fun(ku*gp3[2]+lu,kv*gp3[2]+lv,&x22,&y22,&z22);
            rx+=gc3[0]*gc3[0]*x00+gc3[0]*gc3[1]*x01+gc3[0]*gc3[2]*x02+
                gc3[1]*gc3[0]*x10+gc3[1]*gc3[1]*x11+gc3[1]*gc3[2]*x12+
                gc3[2]*gc3[0]*x20+gc3[2]*gc3[1]*x21+gc3[2]*gc3[2]*x22;
            ry+=gc3[0]*gc3[0]*y00+gc3[0]*gc3[1]*y01+gc3[0]*gc3[2]*y02+
                gc3[1]*gc3[0]*y10+gc3[1]*gc3[1]*y11+gc3[1]*gc3[2]*y12+
                gc3[2]*gc3[0]*y20+gc3[2]*gc3[1]*y21+gc3[2]*gc3[2]*y22;
            rz+=gc3[0]*gc3[0]*z00+gc3[0]*gc3[1]*z01+gc3[0]*gc3[2]*z02+
                gc3[1]*gc3[0]*z10+gc3[1]*gc3[1]*z11+gc3[1]*gc3[2]*z12+
                gc3[2]*gc3[0]*z20+gc3[2]*gc3[1]*z21+gc3[2]*gc3[2]*z22;
        }
    }
    return (rx+ry+rz)*ku*kv;
}

//以下为验证函数
double test(double x)
{
    return x*x*x*x;
}

void test2(double u,double v,double *x,double *y,double *z)
{
    *x=1/(2.0*sqrt(u+v));
    *y=0.0;
    *z=0.0;
}

void test3(double u,double v,double *x,double *y,double *z)
{
    *x=(u+v)*(u+v);
    *y=0.0;
    *z=0.0;
}

int main()
{
    printf("%f/n",Guass_Integral(0.0,1.0,test,15));
    printf("%f/n",Guass_Integral(0.0,1.0,0.0,1.0,test2,10,10)-4.0/3.0*(sqrt(2)-1));
    printf("%f/n",Guass_Integral(0.0,1.0,0.0,1.0,test3,10,10));
    return 0;
}

相关文章推荐

Delphi7高级应用开发随书源码

  • 2003年04月30日 00:00
  • 676KB
  • 下载

Delphi7高级应用开发随书源码

  • 2003年04月30日 00:00
  • 676KB
  • 下载

有限元笔记1-高斯积分

欢迎使用Markdown编辑器写博客在有限元计算中,积分运算是非常重要的一种运算,比如刚度矩阵的运算(以二维问题为例): [K]=∫∫A[B]T[D][B]hdxdy \left[ K \right...

高斯积分MATLAB程序

  • 2011年03月31日 14:04
  • 11KB
  • 下载

高等数学:第十章 曲线积分与曲面积分(3)高斯共识、通量、散度、斯托克斯共识、环流量、旋度

§10.6  高斯公式  通量与散度 一、高斯公式 格林公式表达了平面闭区域上的二重积分与其边界曲线上的曲线积分之间的关系,而高斯公式表达了空间闭区域上的三重积分与其边界曲面上的曲面积分之间的关系...

高斯积分算法源代码

  • 2012年12月28日 00:24
  • 522B
  • 下载

高斯一维积分 fortran

  • 2009年11月15日 15:25
  • 301KB
  • 下载

混合高斯背景建模程序分析

(转)混合高斯背景建模程序分析   标签: 杂谈 分类: 图像处理 将高斯建模改成了用一个亮度分量信息建立,但是发现,修改那个权值,还有那...

高斯积分的matlab实现

  • 2014年01月15日 10:21
  • 7KB
  • 下载
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:三点高斯积分程序
举报原因:
原因补充:

(最多只允许输入30个字)