关闭

三点高斯积分程序

标签: funmath.hini
3476人阅读 评论(0) 收藏 举报
分类:
//高斯积分程序,高斯积分具有计算速度快,精度高,能计算反常积分等优点
#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;
}
0
0

查看评论
* 以上用户言论只代表其个人观点,不代表CSDN网站的观点或立场
    个人资料
    • 访问:333103次
    • 积分:3979
    • 等级:
    • 排名:第8178名
    • 原创:64篇
    • 转载:67篇
    • 译文:2篇
    • 评论:114条
    最新评论
    Opengl网站
    其它
    优秀编程网址