在下面的這段代碼中，包含了高斯-勒讓德、高斯-切比雪夫、以及拉蓋爾和埃爾米特型求積公式，它們分別對應了不同的被積積分型
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