数值计算&Matlab——高斯型积分
Gauss求积
由概念出发，取n个插值节点时，Gauss公式具有2*n-1次代数精度，即对
只看概念会不会很迷，没错，木木刚开始学的时候也是感到很抽象，于是寻找有趣的学习方法使抽象的内容具象化，达到拿起来就可以用的程度，毕竟不是数学专业的学生还需要懂得共识背后的意义，作为一个合格的工科生应该学会用一些有效的数学方法，于是在网上、书上、论坛上搜集各种代码进行分析，使用，废话不多说，上代码！
Matlab程序如下：
function varargout=gausscoef(a,b,n)% 高斯求积公式系数求解if nargin<3n=2; % 默认高斯点endx=sym('x%d',[1,n],'real'); % 节点A=sym('A%d',[1,n],'real'); % 系数eq=sym(zeros(2*n,1)); % 存储方程组for k=1:2*neq(k,1)=dot(A,x.^(k-1)-(b^k-a^k)/k); % 构造方程endvars=symvar(eq); % 提取eq的符号变量fun=matlabFunction(eq,'vars',{vars}); % 将eq转换成匿名函数的形式options=optimoptions(@fsolve,'FunctionTolerance',1e-12);x0=fsolve(fun,2*rand(1,2*n)-1,options); % 求Gauss点和Gauss系数Ac=x0(1:n); % Gauss系数x0(1:n)=[]; % Gauss点[varargout{1:2}]=deal(x0,... % 第一个输出参数为Gauss点Ac); % 第二个输出参数为Gauss求积系数
例1：
a = 2;b = 5;n=2;[x A]=gausscoef(a,b,n)
常见Gauss求积公式
在上式中如果重复计算多次，你会发现得出的结果也会不尽相同，但同样都是高斯点和高斯系数。对于一些特定的形式，一批优秀的科学家给出了相应的精确解析方法，下面以最常见的Gauss-Legendre求积公式（有限元中经常使用）为例。
Gauss-Chebyshev求积、Gauss-Laguerre求积、Gauss-Hermite求积、Gauss-Lobatto求积正交多项式的递推公式以及求积系数，就不在这里一一列举了，教材中大多会提及，老规矩，上代码！
求解正交多项式
function p=opfun(polytype,n,xi)% opfun常用正交多项式% 存储各阶正交多项p=cell(n+1,1);% 递推多项式第二项的系数P1=[1,0;1,0;2,0;-1,1;2,0];% 判断所要求的多项式类型ind=strcmpi(polytype,{'P','T','U','L','H'});% 递推多项式的首项p{1}=1;% 递推多项式的第二项p{2}=P1(ind,:);for k=2:nswitch polytypecase {'P','p'} % 勒让德多项式p{k+1}=(2*k-1)/k*[p{k},0]-(k-1)/k*[0,0,p{k-1}]; % 递推公式case {'T','t','U','u'} % 切比雪夫多项式p{k+1}=2*[p{k},0]-[0,0,p{k-1}]; % 递推公式case {'L','l'} % 拉盖尔多项式p{k+1}=conv([-1,2*k-1],p{k})-(k-1)^2*[0,0,p{k-1}]; % 递推公式case {'H','h'} % 埃尔米特多项式p{k+1}=2*[p{k},0]-2*(k+1)*[0,0,p{k-1}]; % 递推公式otherwise % 如果多项式类型不是以上5种情形，则给出错误提示error('无效的正交多项式类型.')endendif nargin==3 % 如果输入参数有3个for k=1:n+1p_val(k,:)=polyval(p{k},xi); % 根据求得的正交多项式计算指定点处的值endp=p_val; % 将p_val赋给变量pend
求解高斯点和高斯系数
function varargout = Gauss_1(n,type)% 常用高斯求积公式的高斯点及求积系数if nargin<2type = 'P'; % 默认求Gauss-Legendre求积节点和求积系数endp=opfun(type(1),n); % 由type指定正交多项式的系数p=p{end}; % 正交多项式P_n+1(x)x=roots(p); % 求正交多项式的零点，即Gauss点switch typecase {'P',p} % Gauss-Legendre求积A=2./((1-x.^2).*polyval(polyder(p),x).^2); % 求求积系数case {'T','t'} % Gauss-Chebyshev求积A=pi/n*ones(size(x)); % 求求积系数case {'L','l'} % Gauss-Laguerre求积A=gamma(n+1)^2./(x.*polyval(polyder(p),x).^2); % 求求积系数case {'H','h'} % Gauss-Hermite求积A=2^(n+1)*gamma(n+1)*sqrt(pi)./polyval(polyder(p),x).^2;case {'Pl','pl'} % Gauss-Lobatto求积x=roots(polyder(p)); % 求P'(x)的零点A=[2/(n*(n+1));2./(n*(n+1)*polyval(p,x).^2);2./(n*(n+1))];x=[-1;x;1]; % Gauss点otherwiseerror('本函数不能求其他类型的Gauss点及求积系数')end[varargout{1:2}]=deal(x,... % 第一个输出参数为高斯点A); % 第二个输出参数为求积系数
例2：
使用Gauss-Legendre求积公式、Gauss-Chebyshev求积公式，计算积分
f1 = @(x)x.^2./sqrt(1-x.^2);f2 = @(x)x.^2;[x01,A1]=Gauss_1(3,'P');[x02,A2]=Gauss_1(3,'T');I1=dot(A1,f1(x01)) % 高斯-勒让德求积I2=dot(A2,f2(x02)) % 高斯-切比雪夫求积
疑问：在这里大家可能会产生个疑问，你这区间都是[-1,1]，那我实际情况肯定不止于此，那该怎么办呢？
答：记得当时考研的时候，木木看的是张宇的《十八讲》，里面提到积分区间变换的方法，即：
对于一般区间[a,b]上的积分
例3：
利用5点Gauss-Lobatto求积公式，计算积分
a=1;b=5;[x0,A]=Gauss_1(4,'Pl');f3=@(x) sin(x)./x;I3=(b-a)/2*dot(A,f3(x0*(b-a)/2+(a+b)/2)) % Gauss- Lobatto求积
Mathematica
这是一款符号计算特别强大的数学软件，木木正在学习，以后的推送内容也会涉及相关内容，不便于复 制代码，我就直接截图了哈~
