define_gauss_rules.m

matlab实例

function [GaussRule] = define_Gauss_rules
%     SET COORDINATES AND WEIGHTS FOR GAUSSIAN INTEGRATION RULER
%
%     For each rule we have:
%
%     GaussRule(irl).loc = coordinates of integration poins in the interval (-1,1)
%                   .wgt = weights for integration point
%     =================================================================

%     ------------------
%     N = 1
%     ------------------

      GaussRule(1).loc(1) = 0.0e0;

      GaussRule(1).wgt(1) = 2.0e0;

%     ------------------
%     N = 2
%     ------------------

      GaussRule(2).loc(1) = - 0.577350269189626e0;
      GaussRule(2).loc(2) =   0.577350269189626e0;

      GaussRule(2).wgt(1) =   1.0e0;
      GaussRule(2).wgt(2) =   1.0e0;

%     ------------------
%     N  = 3.
%     ------------------

      A      =   0.774596669241483e0;
      H1     =   0.555555555555556e0;
      H2     =   0.888888888888889e0;
      
      GaussRule(3).loc(1) = - A;
      GaussRule(3).loc(2) =   0.0e0;
      GaussRule(3).loc(3) =   A;

      GaussRule(3).wgt(1) =   H1;
      GaussRule(3).wgt(2) =   H2;
      GaussRule(3).wgt(3) =   H1;

%     ------------------
%     N = 4
%     ------------------

      A      =   0.861136311594053e0;
      B      =   0.339981043584856e0;
      
      GaussRule(4).loc(1) = - A;
      GaussRule(4).loc(2) = - B;
      GaussRule(4).loc(3) =   B;
      GaussRule(4).loc(4) =   A;

      GaussRule(4).wgt(1) =   0.347854845137454e0;
      GaussRule(4).wgt(2) =   0.652145154862546e0;
      GaussRule(4).wgt(3) =   GaussRule(4).wgt(2);
      GaussRule(4).wgt(4) =   GaussRule(4).wgt(1);
      
%     ------------------
%     N = 5
%     ------------------

      GaussRule(5).loc(1) = - 0.906179845938664e0;
      GaussRule(5).loc(2) = - 0.538469310105683e0;
      GaussRule(5).loc(3) =   0.0e0;
      GaussRule(5).loc(4) = - GaussRule(5).loc(2);
      GaussRule(5).loc(5) = - GaussRule(5).loc(1);

      GaussRule(5).wgt(1) =  0.236926885056189e0;
      GaussRule(5).wgt(2) =  0.478628670499366e0;
      GaussRule(5).wgt(3) =  0.568888888888889e0;
      GaussRule(5).wgt(4) =  GaussRule(5).wgt(2);
      GaussRule(5).wgt(5) =  GaussRule(5).wgt(1);


      % checking the rules:
      for irl=1:5
          sumLoc = sum( GaussRule(irl).loc );
          sumWht = sum( GaussRule(irl).wgt );
          
          if abs(sumLoc) > 1.0e-14
              error('define_Gauss_rules: wrong rule, loc are not valid');
          end
          if abs(sumWht-2.0) > 1.0e-14
              error('define_Gauss_rules: wrong rule, wht are not valid');
          end
      end