在下面的这段代码中，包含了高斯-勒让德、高斯-切比雪夫、以及拉盖尔和埃尔米特型求积公式，它们分别对应了不同的被积积分型
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\b\'; 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 = X\B\'; 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 = X\B\'; 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 = X\B\'; 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