Numerical integration of $\exp(-x^2)$ in a bounded region

128 Views Asked by At

Is there a possible way to integrate $\exp(-x^2)$ in a bounded 2D triangular region numerically with minimal number of Gauss points? The Gauss-Hermite quadrature scheme is suitable for unbounded regions. I have tried Gauss-Legendre scheme but with 3~5 points the results is not accurate.

2

There are 2 best solutions below

3
On BEST ANSWER

You can use high order Gaussian quadrature rules for triangles. Order $6$ and $7$ quadratures follow (quadrature points and weights correspond to the reference triangle ($[0,0], [1,0], [0,1]$, on other triangles linear transformation has to be used).

if (order == 6)
{
    Allocate(12);
    r[++i] = 0.24928674517091; s[i] = 0.24928674517091; weight[i] = 0.11678627572638 * 0.5;
    r[++i] = 0.24928674517091; s[i] = 0.50142650965818; weight[i] = 0.11678627572638 * 0.5;
    r[++i] = 0.50142650965818; s[i] = 0.24928674517091; weight[i] = 0.11678627572638 * 0.5;
    r[++i] = 0.06308901449150; s[i] = 0.06308901449150; weight[i] = 0.05084490637021 * 0.5;
    r[++i] = 0.06308901449150; s[i] = 0.87382197101700; weight[i] = 0.05084490637021 * 0.5;
    r[++i] = 0.87382197101700; s[i] = 0.06308901449150; weight[i] = 0.05084490637021 * 0.5;
    r[++i] = 0.31035245103378; s[i] = 0.63650249912140; weight[i] = 0.08285107561837 * 0.5;
    r[++i] = 0.63650249912140; s[i] = 0.05314504984482; weight[i] = 0.08285107561837 * 0.5;
    r[++i] = 0.05314504984482; s[i] = 0.31035245103378; weight[i] = 0.08285107561837 * 0.5;
    r[++i] = 0.63650249912140; s[i] = 0.31035245103378; weight[i] = 0.08285107561837 * 0.5;
    r[++i] = 0.31035245103378; s[i] = 0.05314504984482; weight[i] = 0.08285107561837 * 0.5;
    r[++i] = 0.05314504984482; s[i] = 0.63650249912140; weight[i] = 0.08285107561837 * 0.5;
}
else if (order == 7)
{
    Allocate(13);
    r[++i] = 0.33333333333333; s[i] = 0.33333333333333; weight[i] = -0.14957004446768 * 0.5;
    r[++i] = 0.26034596607904; s[i] = 0.26034596607904; weight[i] = 0.17561525743321 * 0.5;
    r[++i] = 0.26034596607904; s[i] = 0.47930806784192; weight[i] = 0.17561525743321 * 0.5;
    r[++i] = 0.47930806784192; s[i] = 0.26034596607904; weight[i] = 0.17561525743321 * 0.5;
    r[++i] = 0.06513010290222; s[i] = 0.06513010290222; weight[i] = 0.05334723560884 * 0.5;
    r[++i] = 0.06513010290222; s[i] = 0.86973979419557; weight[i] = 0.05334723560884 * 0.5;
    r[++i] = 0.86973979419557; s[i] = 0.06513010290222; weight[i] = 0.05334723560884 * 0.5;
    r[++i] = 0.31286549600487; s[i] = 0.63844418856981; weight[i] = 0.07711376089026 * 0.5;
    r[++i] = 0.63844418856981; s[i] = 0.04869031542532; weight[i] = 0.07711376089026 * 0.5;
    r[++i] = 0.04869031542532; s[i] = 0.31286549600487; weight[i] = 0.07711376089026 * 0.5;
    r[++i] = 0.63844418856981; s[i] = 0.31286549600487; weight[i] = 0.07711376089026 * 0.5;
    r[++i] = 0.31286549600487; s[i] = 0.04869031542532; weight[i] = 0.07711376089026 * 0.5;
    r[++i] = 0.04869031542532; s[i] = 0.63844418856981; weight[i] = 0.07711376089026 * 0.5;
}
1
On

A quadrature rule for tetrahedrons of order 5 (taken from enter link description here, page 12):

if (order == 5) //15-point quadrature
{
    r = new double[15];
    s = new double[15];
    t = new double[15];
    weight = new double[15];
                    
    double[] w = new double[] {0.030283678097089, 0.006026785714286, 0.011645249086029, 0.010949141561386};
    double[] l1 = new double[] {0.0, 0.333333333333333333};
    double[] l2 = new double[] {0.727272727272727272, 0.090909090909090909};
    double[] l3 = new double[] {0.066550153573664, 0.433449846426336};
    
    int i = -1;

    i++; r[i] = 0.25; s[i] = 0.25; t[i] = 0.25; weight[i] = w[0];
    
    i++; r[i] = l1[1]; s[i] = l1[1]; t[i] = l1[1]; weight[i] = w[1];
    i++; r[i] = l1[0]; s[i] = l1[1]; t[i] = l1[1]; weight[i] = w[1];
    i++; r[i] = l1[1]; s[i] = l1[0]; t[i] = l1[1]; weight[i] = w[1];
    i++; r[i] = l1[1]; s[i] = l1[1]; t[i] = l1[0]; weight[i] = w[1];
    
    i++; r[i] = l2[1]; s[i] = l2[1]; t[i] = l2[1]; weight[i] = w[2];
    i++; r[i] = l2[0]; s[i] = l2[1]; t[i] = l2[1]; weight[i] = w[2];
    i++; r[i] = l2[1]; s[i] = l2[0]; t[i] = l2[1]; weight[i] = w[2];
    i++; r[i] = l2[1]; s[i] = l2[1]; t[i] = l2[0]; weight[i] = w[2];
    
    i++; r[i] = l3[0]; s[i] = l3[1]; t[i] = l3[1]; weight[i] = w[3];
    i++; r[i] = l3[1]; s[i] = l3[0]; t[i] = l3[1]; weight[i] = w[3];
    i++; r[i] = l3[1]; s[i] = l3[1]; t[i] = l3[0]; weight[i] = w[3];
    i++; r[i] = l3[1]; s[i] = l3[0]; t[i] = l3[0]; weight[i] = w[3];
    i++; r[i] = l3[0]; s[i] = l3[1]; t[i] = l3[0]; weight[i] = w[3];
    i++; r[i] = l3[0]; s[i] = l3[0]; t[i] = l3[1]; weight[i] = w[3];
}