Gaussian Quadrature rule for Triangle and Tetrahedron Qikun Wu, Liuxing Shen Introduction We demonstrate the position of the Gaussian points in 2D and 3D case, and finished task 2.

1. two dimension case in two dimension case, we will talk about square and triangle. 1.1 tensor products of the one dimensional formula intuitively, if we directly use tensor product of one-dimensional case, we can have a ‘Gaussian Quadrature’ scheme in two dimension, which is illustrated as follows:

 

1

1

1 1 1

f ( , )d d  

1

1



1

1

f ( , )d  d



n 1  n      H i f (i , ) d   H i  f (i , )d 1 1  i 1  i 1 n  n  n n   H i   H j f (i , j )    H i H j f (i , j ) i 1  j 1  i 1 j 1

And this tensor product can extend to arbitrary dimension:

  n 1

1

1 1 n



1

1

f ( x1 , x2 n i1

xm )d x1dx2

dxm xm( im ) )

  i1 1 i2 1

im 1

H

Hi2

Him f ( x1( i1 ) , x2( i2 )

However, we can see that this quadrature is not based on some ‘orthogonal polynomial in 2D’, which is not an optimal solution, so we can expect that we can choose less point to reach the same accuracy. So we should let the quadrature scheme with the following form:

I 

1

1 1



1

f ( , )d d  Wi f (i ,i ) i 1

n

1.2 non-tensor-product formula As stated before, the non-tensor product form can be made using only 7 points to achieve the same accuracy while the tensor form need 9 points. However, since the orthogonal polynomials are unknown in two and three dimension, these non-tensor-product form are complicated to derive. The corresponding coefficients are usually determined by method of underdetermined coefficients. We shall show the example of triangle in the next section. 1.3 quadrature rule in triangle We now discuss the scenario where we will do integral on a canonical triangle, where the coordinates of the vertexes are [0,0], [0,1], [1,0]. We first consider the simplest case:

 f ( , )d d  W f ( , )  E
1 1 1 

We shall expect that the quadrature rule is exact for the monomials 1,  , , which means:


0

1

1

0 1

1d  d  W1

1  W11 0 0 6 1 1 1 0 0 d d  6  W11



1

 d  d 

The solution to the corresponding formulas is explicit:

1 1 W1  , 1  1  2 3
So the one-point quadrature rule is:

 f ( , )d d  2 f ( 3 , 3)  E


1

1 1

In which we can see that ( , ) is exactly the centroid of the triangle. For an example with three points, we also give the formula:

1 1 3 3

 f ( , )d d  6  f ( 6 , 6 )  f ( 6 , 6 )  f ( 6 , 6 )  E  


1

1 1

4 1

1 4 

For an example with four points, the formula is:

 f ( , )d d   96 f ( 3 , 3)  96  f ( 5 , 5 )  f ( 5 , 5 )  f ( 5 , 5 )   E  


27

1 1

25 

1 1

3 1

1 3 

For a general case, the quadrature rule is expressed as:




i f ( , )d d  Ae Wi f ( 1i , 2 , 3i )  E i 1

n

Where Ae is the area of the triangle and ( 1 ,  2 ,  3 ) is the triangular coordinates, which i i i

means the position of each point can be represented as the area constructed by the point and one edge divided by the area of the triangle. For instance, if the coordinate is

1 1 1  , ,  , then this point is exactly the centroid of the triangle. In the canonical triangle,  3 3 3 we can simply set

 2i  i , 3i  i .

The following table, given by D.A. Dunavant1, lists the triangle coordinates of different quadrature rule, which can approximate the polynomial up to 20.
1

D.A. Dunavant, 1985, High Degree Efficient Symmetrical Gaussian quadrature rules for the triangle. International Journal of Numerical Methods in Engineering, vol. 21, pp.1129-1148

2. three dimension case From the general expression of the quadrature rule in two dimensional case, we can write the general expression of the corresponding rule in three dimensional case:

 f ( x, y, z )dxdydz  V W f ( e  i 1 i

n

i 1

i i ,  2 ,  3i ,  4 )  E

Where for three dimensional regions, ( 1 ,  2 ,  3 ,  4 ) represents the volume coordinate, which i i i i

means the position of each point can be represented as the volume constructed by the point and one of the surface divided by the total volume Ve . The quadrature rules to order six and order eight are introduced by Jinyun2 and Keast3. Similarly, in the case that the vertexes of the tetrahedron is (0,0,0), (1,0,0), (0,1,0), (0,0,1), we just let

 2i  xi , 3i  yi , 4i  zi .

3. illustration (task 1) we plot the point according to the table by D.A. Dunavant4 and Keast5. The gp2d.m and gp3d.m is the corresponding MATLAB code

2

Y. JInyun, 1984, Symmetric Gaussian quadrature formulae for tetrahedronal regions. Computer Methods in Applied Mechanics an Engineering, 43:349-353. 3 P. Keast, 1984, Moderate-degree tetrahedral quadrature formulas. Computer methods in Applied Mechanics and Engineering, 55:339-348 4 D.A. Dunavant, 1985, High Degree Efficient Symmetrical Gaussian quadrature rules for the triangle. International Journal of Numerical Methods in Engineering, vol. 21, pp.1129-1148 5 P. Keast, 1984, Moderate-degree tetrahedral quadrature formulas. Computer methods in Applied Mechanics and Engineering, 55:339-348

4. Task 2 We now deal with the integral

If  

1

0

e
0

1

x y

sin( x  y )dxdy

First, MATLAB function quad2d gives the integration value of 0.843808868002693, with estimated error boundary of 6.980097461608347e-11. So below, we can view the integration value as the exactly value. We get the numerical integration value by implementing Gaussian point in 1D twice, and the results as below. For more detail about the method, check out code file g2d.m.

Gaussian point 1 4 9 16

numerical integration 0.841470985 0.845366941 0.843809323 0.843808853

error 0.002337883 0.001558073 4.55E-07 1.51E-08

Interesting results can be seen in the table, this method has the property that second order convergence alternates with first order convergence, which is equivalent to 2 order convergence in general.

Appendices

gp2d.m (Task 1, 2D case)
% generate visualization of gaussian point in 2D x=[0;0;1]; y=[1;0;0]; fill(x,y,'w','MarkerSize',20); hold on; % p=1 ng=1 % x1=1/3; y1=1/3; % plot(x1,y1,'r.','MarkerSize',30); % title('plot of Gaussian point with precision 1'); % p=2 ng=3 % x2=[1/6,1/6,2/3]; % y2=[1/6,2/3,1/6]; % plot(x2,y2,'b.','MarkerSize',20); % title('plot of Gaussian point with precision 2'); % p=3 ng=4 % x3=[0.6,0.2,0.2]; % y3=[0.2,0.6,0.2]; % plot(x3,y3,'g.','MarkerSize',20); % plot(1/3,1/3,'g.','MarkerSize',20); % title('plot of Gaussian point with precision 3'); % hold off;

% p=4 ng=6 % x4=[0.108103018168070,0.445948490915965,0.445948490915965]; % y4=[0.445948490915965,0.108103018168070,0.445948490915965]; % plot(x4,y4,'b.','MarkerSize',20); % x4=[0.816847572980459,0.091576213509771,0.091576213509771]; % y4=[0.091576213509771,0.091576213509771,0.816847572980459]; % plot(x4,y4,'b.','MarkerSize',20); % title('plot of Gaussian point with precision 4'); % p=5 ng=7 % plot(1/3,1/3,'b.','MarkerSize',20); % x5=[0.059715871789770,0.470142064105115,0.470142064105115]; % y5=[0.470142064105115,0.059715871789770,0.470142064105115];

%

plot(x5,y5,'b.','MarkerSize',20);

% x5=[0.797426985353087,0.101286507323456,0.101286507323456]; % y5=[0.101286507323456,0.797426985353087,0.101286507323456]; % plot(x5,y5,'b.','MarkerSize',20); % title('plot of Gaussian point with precision 5');

% p=6 ng=12 x6=[0.501426509658179,0.249286745170910,0.249286745170910]; y6=[0.249286745170910,0.501426509658179,0.249286745170910]; plot(x6,y6,'m.','MarkerSize',20); x6=[0.873821971016996,0.063089014491502,0.063089014491502]; y6=[0.063089014491502,0.873821971016996,0.063089014491502]; plot(x6,y6,'m.','MarkerSize',20); x6=[0.053145049844817,0.053145049844817,0.310352451033784,0.310352451 033784,0.636502499121399,0.636502499121399]; y6=[0.310352451033784,0.636502499121399,0.053145049844817,0.636502499 121399,0.053145049844817,0.310352451033784]; plot(x6,y6,'m.','MarkerSize',20); title('plot of Gaussian point with precision 6');

gp3d.m (Task 1, 3D case)
% generate visualization of gaussian point in 3D x=[0,1,0,0,0,0]; y=[0,0,1,0,0,1]; z=[0,0,0,0,1,0]; plot3(x,y,z,'MarkerSize',20); hold on; x=[0,1]; y=[0,0]; z=[1,0]; plot3(x,y,z,'MarkerSize',20); %x=[1,0,0,1,1,1,0,0,1,1]; %y=[0,0,1,1,0,0,0,1,1,0]; %z=[1,1,1,1,1,0,0,0,0,0]; %plot3(x,y,z,'MarkerSize',20); %x=[1,1]; %y=[1,1]; %z=[1,0]; %plot3(x,y,z,'MarkerSize',20); % n=1 p=1 % x=1/4;

% y=1/4; % z=1/4; % scatter3(x,y,z,200,'black.'); % title('plot of Gaussian point with precision 1');

% n=4 p=2 % x=[0.585410196624969,0.138196601125011,0.138196601125011,0.1381966011 25011]; % y=[0.138196601125011,0.585410196624969,0.138196601125011,0.1381966011 25011]; % z=[0.138196601125011,0.138196601125011,0.585410196624969,0.1381966011 25011]; % scatter3(x,y,z,200,'b.'); % title('plot of Gaussian point with precision 2'); % n=5 p=3 % x=1/4; % y=1/4; % z=1/4; % scatter3(x,y,z,400,'g.'); % x=[0.5,1/6,1/6,1/6]; % y=[1/6,0.5,1/6,1/6]; % z=[1/6,1/6,0.5,1/6]; % scatter3(x,y,z,200,'g.'); % title('plot of Gaussian point with precision 3'); % % n=11 p=4 x=1/4; y=1/4; z=1/4; scatter3(x,y,z,800,'m.'); x=[0.785714285714286,0.071428571428571,0.071428571428571,0.0714285714 28571]; y=[0.071428571428571,0.785714285714286,0.071428571428571,0.0714285714 28571]; z=[0.071428571428571,0.071428571428571,0.785714285714286,0.0714285714 28571]; scatter3(x,y,z,200,'m.'); x=[0.399403576166799,0.399403576166799,0.399403576166799,0.1005964238 33201,0.100596423833201,0.100596423833201];

y=[0.399403576166799,0.100596423833201,0.100596423833201,0.3994035761 66799,0.100596423833201,0.399403576166799]; z=[0.100596423833201,0.399403576166799,0.100596423833201,0.1005964238 33201,0.399403576166799,0.399403576166799]; scatter3(x,y,z,200,'m.'); title('plot of Gaussian point with precision 4');

g2d.m (Task 2)
% two linear gauss numerical integration to approximate 2d integration gm=zeros(5,5); gm(2,1)=-sqrt(3)/3; gm(2,2)=sqrt(3)/3; gm(3,1)=-sqrt(3/5); gm(3,3)=sqrt(3/5); gm(4,1)=-sqrt((3+2*sqrt(6/5))/7); gm(4,2)=-sqrt((3-2*sqrt(6/5))/7); gm(4,3)=sqrt((3-2*sqrt(6/5))/7); gm(4,4)=sqrt((3+2*sqrt(6/5))/7); gm(5,1)=-(1/3)*sqrt(5+2*sqrt(10/7)); gm(5,2)=-(1/3)*sqrt(5-2*sqrt(10/7)); gm(5,4)=(1/3)*sqrt(5-2*sqrt(10/7)); gm(5,5)=(1/3)*sqrt(5+2*sqrt(10/7)); re=zeros(4,1); syms x y f(x,y); f(x,y)=exp(x-y)*sin(x+y); % one point re(1,1)=eval(f(1/2,1/2)); % two points y1=(gm(2,1)-(-1))/2; y2=(gm(2,2)-(-1))/2; re(2,1)=1/4*(eval(f(y1,y1))+eval(f(y1,y2))+eval(f(y2,y1))+eval(f(y2,y 2))); % three points y1=(gm(3,1)-(-1))/2; y2=(gm(3,2)-(-1))/2; y3=(gm(3,3)-(-1))/2; re(3,1)=5/18*(5/18*eval(f(y1,y1))+8/18*eval(f(y2,y1))+5/18*eval(f(y3, y1))); re(3,1)=re(3,1)+8/18*(5/18*eval(f(y1,y2))+8/18*eval(f(y2,y2))+5/18*ev

al(f(y3,y2))); re(3,1)=re(3,1)+5/18*(5/18*eval(f(y1,y3))+8/18*eval(f(y2,y3))+5/18*ev al(f(y3,y3))); % four points y1=(gm(4,1)-(-1))/2; y2=(gm(4,2)-(-1))/2; y3=(gm(4,3)-(-1))/2; y4=(gm(4,4)-(-1))/2; a=(18-sqrt(30))/72; b=(18+sqrt(30))/72; re(4,1)=a*(a*eval(f(y1,y1))+b*eval(f(y2,y1))+b*eval(f(y3,y1))+a*eval( f(y4,y1))); re(4,1)=re(4,1)+b*(a*eval(f(y1,y2))+b*eval(f(y2,y2))+b*eval(f(y3,y2)) +a*eval(f(y4,y2))); re(4,1)=re(4,1)+b*(a*eval(f(y1,y3))+b*eval(f(y2,y3))+b*eval(f(y3,y3)) +a*eval(f(y4,y3))); re(4,1)=re(4,1)+a*(a*eval(f(y1,y4))+b*eval(f(y2,y4))+b*eval(f(y3,y4)) +a*eval(f(y4,y4))); syms x y ; f=@(x,y) exp(x-y).*sin(x+y); [I,errbnd]=quad2d(f,0,1,0,1);