数值积分——高斯型求积公式

　　在下面的这段代码中，包含了高斯-勒让德、高斯-切比雪夫、以及拉盖尔和埃尔米特型求积公式，它们分别对应了不同的被积积分型

　　1.代码

%%高斯型求积公式
%%Y是函数表达式，interval是求积区间，n是求积阶数
%%对于求一般形式的非反常积分，可用勒让德型，
%%对于求形如f(x)/sqrt(1-x^2)的非反常积分，可用第一类切比雪夫型，
%对于形如f(x)*sqrt(1-x^2)的非反常积分，可用第二类切比雪夫型，切比雪夫型积分应在[-1 1]上
%对于反常积分f(x)*exp(-x)或者f(x)*exp(-x^2)，且区间为[0,+inf]或[-inf,+inf]，可用拉盖尔或者埃尔米特型（注意Y是f(x)）
function GQF = Gauss_Quadrature_formula(Y,interval,n,type)
x = sym('x');
global sum;

%%勒让德型
if strcmp(type,'Legendre') == 1
    a = interval(1);b = interval(2);
    s = (b-a)*x/2+(a+b)/2;
    F = subs(Y,x,s);
    X0 = solve([Orthogonal_polynomial(type,n+1)],[0])';
    for k = 1:n
        X(k,:) = X0.^(k-1);
        B(k) = int(x^(k-1),-1,1);
    end
    A = X(B');
    F_value = subs(F,X0);
    sum = 0;
    for k = 1:n
        sum = sum+A(k)*F_value(k);
    end
    sum = (b-a)*sum/2;
    
    %%拉盖尔型
elseif strcmp(type,'Laguerre') == 1
    X0 = solve([Orthogonal_polynomial(type,n+1)],[0])';
    for i = 1:n
        X(i,:) = X0.^(i-1);
        b(i) = factorial(i-1);
    end
    A = X';
    F_value = subs(Y,X0);
    sum = 0;
    for i = 1:n
        sum = sum+A(i)*F_value(i);
    end
    
    %%埃尔米特型
elseif strcmp(type,'Hermite') == 1
    X0 = solve([Orthogonal_polynomial(type,n+1)],[0])';
    for k = 1:n
        X(k,:) = X0.^(k-1);
        if (ceil(k/2) == k/2) == 1
            B(k) =0;
        else
            B(k) = gamma(k/2);
        end
    end
    A = XB';
    F_value = subs(Y,X0);
    sum = 0;
    for k = 1:n
        sum = sum+A(k)*F_value(k);
    end
    
    %%切比雪夫型
elseif strcmp(type(1:9),'Chebyshev') == 1
    class = eval(type(10));
    if class == 1
        X0 = solve([Orthogonal_polynomial(type,n+1)],[0])';
        for k = 1:n
            X(k,:) = X0.^(k-1);
            if (ceil(k/2) == k/2) == 1
                B(k) =0;
            else
                h = k-1;
                B(k) = pi*Double_factorial(h+1)/(Double_factorial(h)*(h+1));
            end
        end
        A = XB';
        F_value = subs(Y,X0);
        sum = 0;
        for k = 1:n
            sum = sum+A(k)*F_value(k);
        end
    elseif class == 2
        X0 = solve([Orthogonal_polynomial(type,n+1)],[0])';
        for k = 1:n
            X(k,:) = X0.^(k-1);
            if (ceil(k/2) == k/2) == 1
                B(k) =0;
            else
                h = k-1;
                B(k) = pi*Double_factorial(h+1)/(Double_factorial(h+2)*(h+1));
            end
        end
        A = XB';
        F_value = subs(Y,X0);
        sum = 0;
        for k = 1:n
            sum = sum+A(k)*F_value(k);
        end
    end
end
GQF = vpa(sum,max([2,2^min([n,3])]));

%%组合数中规定n>=m
    function NC = Number_of_combinations(m,n)
        NC = factorial(n)/(factorial(m)*factorial(n-m));
        function F = factorial(n)
            if n == 0
                F = 1;
            else
                F = factorial(n-1)*n;
            end
        end
    end

%%双阶乘
    function DF = Double_factorial(n)
        if n == 0
            DF = 1;
        elseif n == 1
            DF = 1;
        elseif n == -1
            DF = -1;
        elseif n > 1
            DF = Double_factorial(n-2)*n;
        elseif n < -1
            DF = Double_factorial(n+2)*n;
        end
    end

%%阶乘函数
    function F = factorial(n)
        if n == 0
            F = 1;
        else
            F = factorial(n-1)*n;
        end
    end

end

　　

%%正交多项式
%此函数包括勒让德正交多项式（定义区间[-1,1]），切比雪夫正交多项式（两类，
%在这里，规定第一类切比雪夫多项式是以1/sqrt(1-x^2)作为权函数，第二类切比雪夫多项式以sqrt(1-x^2)作为权函数得到的）（定义区间[-1,1]）
%拉盖尔正交多项式（定义区间[0 +inf]），埃尔米特正交多项式（定义区间[-inf +inf]），输入项数N应从1开始
%%n是多项式的项数，n>=0,type是类型，分为Legendre、Chebyshev1、Chebyshev2、Laguerre、Hermite以及幂函数{1,x,x^2,x^3,…}对应其正交多项式
function OP = Orthogonal_polynomial(type,N)
sym type;
if strcmp(type,'Legendre') == 1
    L = Legendre(N);
    OP = simplify(L(N));
elseif strcmp(type,'Hermite') == 1
    H = Hermite(N);
    OP = simplify(H(N));
elseif strcmp(type,'Laguerre') == 1
    La = Laguerre(N);
    OP = simplify(La(N));
elseif strcmp(type,'幂函数') == 1
    Po = Power_fun(N)
    OP = simplify(Po(N));
elseif strcmp(type(1:9),'Chebyshev') == 1
    class = eval(type(10));
    Che = Chebyshve(N,class);
    OP = simplify(Che(N));
end
%%幂函数正交
    function Po = Power_fun(N)
        x = sym('x');
        for i = 1:N
            Power(i) = x^(i-1);
        end
        Po = Power;
    end
%%勒让德多项式
    function L = Legendre(N)
        x = sym('x');
        for i = 1:N
            Leg(i) = diff((x^2-1)^(i-1),i-1)/(factorial(i-1)*2^(i-1));
        end
        L = Leg;
    end
%%切比雪夫多项式
    function C = Chebyshve(n,class)
        x = sym('x');
        if class == 1
            T = string([1 x]);
            T = sym(T);
            if n <=2
                C = T(1:n);
            else
                for i = 2:n
                    T(i+1) = 2*x*T(i)-T(i-1);
                end
                C = T(1:n);
            end
        elseif class ==2
            U = string([1]);
            U = sym(U);
            U = [U 2*x];
            if n <=2
                C = U(1:n);
            else
                for i = 2:n
                    U(i+1) = 2*x*U(i)-U(i-1);
                end
                C = U(1:n);
            end
        end
    end
%%埃尔米特多项式
    function H = Hermite(N)
        x = sym('x');
        for i = 1:N
            He(i) = (-1)^N*exp(x^2)*diff(exp(-x^2),(i-1));
        end
        H = simplify(He);
    end
%%拉盖尔多项式
    function La = Laguerre(N)
        x = sym('x');
        for i = 1:N
            Lag(i) = exp(x)*diff(x^(i-1)*exp(-x),(i-1));
        end
        La = simplify(Lag);
    end
%%阶乘函数
    function F = factorial(n)
        if n == 0
            F = 1;
        else
            F = factorial(n-1)*n;
        end
    end
end

　　2.例子

　　（1）高斯-勒让德求积公式

syms x;
n = 5;

disp('高斯-勒让德求积公式');
Y1 = exp(x)*sin(x)+log(x+1);
interval1=[0 pi];
type1 = 'Legendre';
disp('求积公式值：');
Gauss_Quadrature_formula(Y1,interval1,n,type1)
disp('真实值');
vpa(int(Y1,x,interval1),9)

　　

高斯-勒让德求积公式
求积公式值：
ans =
14.814301
真实值
ans =
14.8142899

　　（2）I型切比雪夫

syms x;
n = 5;

disp('高斯-切比雪夫I型');
Y2 = cos(x^4)*exp(-abs(x));
interval2 = [-1 1];
type2 = 'Chebyshev1';
disp('I型切比雪夫');
Gauss_Quadrature_formula(Y2,interval2,n,type2)
disp('真实值');
vpa(int(Y2/sqrt(1-x^2),x,interval2),9)

　　

高斯-切比雪夫I型
I型切比雪夫
ans =
1.6533496
真实值
ans =
1.58844833

　　（3）II型切比雪夫

syms x;
n = 5;

disp('高斯-切比雪夫II型');
Y3 = cos(x^3);
interval3 = [-1 1];
type3 = 'Chebyshev2';
disp('II型切比雪夫');
Gauss_Quadrature_formula(Y3,interval3,n,type3)
disp('真实值');
vpa(int(Y3*sqrt(1-x^2),x,interval3),9)

　　

高斯-切比雪夫II型
II型切比雪夫
ans =
1.5113594
真实值
ans =
1.51150634

　　(4)拉盖尔型

syms x;
n = 5;

Y4 = x*sin(x+cos(x));
interval4 = [0 inf];
type4 = 'Laguerre';
disp('拉盖尔型');
Gauss_Quadrature_formula(Y4,interval4,n,type4)
disp('真实值');
vpa(int(Y4*exp(-x),x,interval4),9)

　　

拉盖尔型
ans =
0.83601799
真实值
ans =
0.823632836

　　（5）埃尔米特型

syms x;
n = 5;

Y5 = cos(x+sin(x)-1);
interval5 = [-inf inf];
type5 = 'Hermite';
disp('埃尔米特型');
Gauss_Quadrature_formula(Y5,interval5,n,type5)
disp('真实值');
vpa(int(Y5*exp(-x^2),x,interval5),9)

　　

埃尔米特型
ans =
0.40264915
真实值
ans =
0.395314636