Matlab数值积分

function [gc,Ak,xk]=guasschebyshev(fun,n)
% 高斯-勒让德数值积分
%
% 参数说明
% fun：积分表达式，可以是函数句柄、inline函数、匿名函数、字符串表达式，但是必须可以接受矢量输入
% a,b：积分上下限只能[-1 1]，故忽略
% n：积分阶数，可以任意正整数，但是不建议设置过大，大不一定能得到更好的精度
% gc：积分结果
% Ak：求积系数
% xk：求积节点，满足gc=sum(Ak.*fun(xk))
%
% 举例说明
% fun=@(x)exp(x).*cos(x); % 必须可以接受矢量输入
% quadl(fun,-1,1) % 调用MATLAB内部积分函数检验
% [gc,Ak,xk]=guasschebyshev(fun,7)
%
% 高斯积分说明
% 1.高斯积分是精度最高的插值型数值积分，具有2n+1阶精度，并且高斯积分总是稳定。一般插值型数值积分精度为至少n阶，且具有高阶不稳定性
% 2.高斯求积节点就是对应n阶正交多项式的根，构建首项系数为1的正交多项式
% 3.高斯求积系数，可以由Lagrange多项式插值系数进行积分得到
% 4.由高斯求积节点为根构成的多项式，与任何不超过n阶的多项式正交
%
% 勒让德正交多项式说明
% 1.区间[-1,1]上关于权rho(x)=1/sqrt(1-x^2)的正交多项式系，满足
%                                   |-   pi/2      (i=j≠0)
% ∫(Ti(x)*Tj(x)/sqrt(1-x^2),-1,1)= |    pi        (i=j=0)
%                                   |-   0         (i≠j)
% 2.切比雪夫正交多项式的通式为：T0=1, Tn=cos(n*arcos(x))  (n=1,2,...)
% 3.Tn(x)=0在[-1,1]有n个零点xk=cos((2k-1)/(2n)*pi) (k=1,2,...n)，
%                  有n+1个最值点pk=cos(k*pi/n)  (k=0,2,...,n)
% 3.关于高斯-切比雪夫的数值积分的求积系数Ak=pi/(n+1)，求积节点xk=cos((2k+1)/(2n+2)*pi)  (k=0,2,...,n)
%
% 计算求积系数
k=0:n;
Ak=ones(1,n+1)*pi/(n+1);
% 计算求积节点
xk=cos((2*k+1)/(2*n+2)*pi);
% 计算权函数
rho=@(x)sqrt(1-x.^2);
% 积分函数有效性检验
fun=fcnchk(fun,'vectorize');
% 带权的积分函数值
fx=fun(xk).*rho(xk);
% 计算积分值
gc=sum(Ak.*fx);