三点高斯积分程序

原创 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;
}

单元节点和积分点有什么区别

单元节点和积分点有什么区别 学过数值积分的应该知道,有限元中的积分点指高斯积分点,因为这些点的收敛性好,精度高。  r]3SE7xCvM(S'Wy2H  t u2xU:S5L...
  • zhuxianjianqi
  • zhuxianjianqi
  • 2013年06月14日 09:54
  • 8622

三点高斯积分程序

//高斯积分程序,高斯积分具有计算速度快,精度高,能计算反常积分等优点#include #include #define PI 3.14159265static const double gp3[3]...
  • Y___Y
  • Y___Y
  • 2007年09月17日 18:04
  • 3951

简单潮流计算——PQ分解法和牛顿-拉夫逊法

  • 2009年03月03日 21:46
  • 21KB
  • 下载

数值积分之Gauss求积法五点公式

//Gauss求积法五点公式 #include #include using namespace std; class quadrature { private:  int i; ...
  • u011378809
  • u011378809
  • 2013年08月08日 18:38
  • 1540

计算高斯积分与插值

print "xxx",np.sqrt(np.pi) integrate.quad(lambda x: np.exp(-x**2), -np.inf,np.inf) 差值即在数据集...
  • liyanbocs
  • liyanbocs
  • 2016年10月07日 23:36
  • 891

数值计算方法与C语言工程函数库

  • 2012年03月20日 10:39
  • 15.55MB
  • 下载

室内定位常用算法概述

一. 室内定位目的和意义 随着数据业务和多媒体业务的快速增加,人们对定位与导航的需求日益增大,尤其在复杂的室内环境,如机场大厅、展厅、仓库、超市、图书馆、地下停车场、矿井等环境中,常常需...
  • fkepgydhbyuan
  • fkepgydhbyuan
  • 2017年01月08日 10:37
  • 3043

等参元的高斯积分详解

等参元的高斯积分详解
  • DaI253
  • DaI253
  • 2015年12月08日 14:54
  • 146

高斯积分

最近看书,总是提到高斯积分。自己也知道有这么个方法,但是一直不是很清楚。所以这次尝试写了一下,的确明白了一些。但是,写的不好,代码没有优化,也没有扩展,只能算到三次积分,且最多能取5点。以后有机会希望...
  • begtostudy
  • begtostudy
  • 2006年11月15日 16:41
  • 2242

求解偏微分方程开源有限元软件deal.II学习--Step 3

求解偏微分方程开源有限元软件deal.II学习--Step 3 Posted on 2016-08-25   |   In computational material science   |   ...
  • zyex1108
  • zyex1108
  • 2016年08月30日 15:03
  • 952
内容举报
返回顶部
收藏助手
不良信息举报
您举报文章:三点高斯积分程序
举报原因:
原因补充:

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